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BUCKLING OF A TRIANGULAR PLATE 
SHEARING FORCES 


By J. M. KLITCHIEFF (Belgrade) 
[Received 13 February 1951] 


SUMMARY 


The dimensions of built-up plate girders are often taken so as to eliminate buckling 
of the web under service conditions. The corresponding spacing of transverse 
stiffeners is usually determined by considering a portion of the web between two 
stiffeners as a simply supported rectangular plate under the action of shearing forces. 

In addition to transverse stiffeners, the web is often reinforced by diagonal 
stiffening ribs. To estimate their effect it is necessary, however, to know the value 
of the shearing force for which a simply supported triangular plate first begins to 
buckle (Fig.). 

In what follows the spacing of transverse stiffeners is assumed equal to the depth 
of the web, so that the problem is that of the buckling of a right-angled isosceles 


triangle. 











Notation 
v deflexion of plate at point (a, y). 
D its flexural rigidity. 


T shearing force per unit length. 


even integers. 
odd integers. 
constants. 
7* D 
32 aT 
[Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 3 (1951)] 
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Details of solution 

LET us consider a right-angled isosceles triangular plate, simply supported 
along its edges, and acted on by shearing forces 7’, uniformly distributed 
along the edges (Fig.) « = a and y = 0, and by normal forces of the same 
magnitude 7’, uniformly distributed along the edge y = x. The differential 
equation of the deflected middle surface then becomes 


a9 09 9 9 9 
co” Gg cru ofw Oo-w 

D+ Nog t+o3) +27 5 = 0. 
of Oy cx" ~§ oy oxcy 


The boundary conditions for a simply supported plate can be expressed 


in the form 


Of-w 
along the edge x = a: w = 0, = 0; 
z ‘ 9 
ox 
9 
Cw 
along the edge y = 0: w= 0. —~=0; 
ey? 
. . 
o-Ww o-w 
along the edge y = x: w = 0, -~+—.= 0 
Caz § Oy? 


These conditions will be satisfied if the solution of the differential equation 
is taken in the form 


L L 


“ “ - MAX . nw 
w > » C. {sin sin — 
an or a a 


— 
m=2,4.. n=1,3... 


} - nae . mary 
a -Sin sin —}. 
a a 


Substituting in the differential equation 


= — a\? - mre. nary . NaxX . mMmr7y 
S S C {D( (m?+- n?)(sin - pe an —_ | 1 
hood ideal wer Ne a a a a 
m=2,4... n=1,3... 


me MTL 7Tt NTL miTry 
1-2 1 mn{cos cos — 7. cos cos : — (). 
a a a a 


Multiplying this equation by 


A or ee , : 
sin / sin - Y dxdy (2 By Mies 9 Te. Bee 
a a 


integrating with respect to x and y, from 0 to a, and observing that 


a a 


. . ») . > . 
Marx. 7a 2am Mae. jr 
cos sind dx ——.., cos sin’ dx — 0, 
, a a 7 m*—)/* a a 
0 0 


there results 


x 
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The condition that these equations should have a non-zero solution gives 


the infinite determinantal equation 


5°A ; —f 35 63 io 
13°A+-35 55 19 83 
B Ay —H 9 
3 ¢ 4 BM. =O. 
‘ 2 9 Me GPA+ He 


If the first row and the first column of this determinant is divided by 5, 
the second row and the second column by 13, and so on, the equation 
becomes a secular equation of symmetrical type, having only real roots, 
the largest of which determines the critical value of the shearing force 7... 

The calculation can be effected by a process of successive approximation. 
The determinantal equation formed by ignoring all elements except those 
common to the first three rows and columns gives A = 0-0047; including 
elements common to the first four, five, and seven rows and columns gives 


\ 00-0061. A 0-00666. and A 0-00666 respectively.+ 


Conclusion 


The magnitude of the critical shearing force corresponding to A ~ 0-0067 


i D 
re 467° 
' a" 
since the value of the critical shearing force applied to a square plate is (1) 
) 
1 Y-4 ry. . 
a> 


it follows that the effect on the value of the critical shearing force of 


reinforcing the web by additional diagonal stiffeners is considerable. 


Cr} ts were c: ed for the writer by Mrs. I. Yassensky. 
| REFERENCE 
1. S. Timoshenko, I Stability (New York, 1933), p. 360. 








ON THE TWO-DIMENSIONAL PROBLEM 
OF A SEMI-INFINITE ELASTIC BODY 
COMPRESSED BY AN ELASTIC PLANE 


By H. OKUBO 


(Institute of High Speed Mechanics, 
Tohoku University, Sendai, Japan) 


[Received 21 March 1950] 


SUMMARY 


An approximate solution of the two-dimensional problem of an elastic foundation 
compressed by a pillar is given, account being taken of the elasticity of the pillar. 
It is found that the contact pressure between the two elastic planes in contact has 


a tendency to concentrate at both edges of the prism. 


1. THE object of this paper is to find an approximate solution of the two- 
dimensional problem of an elastic foundation compressed by a pillar with 
a rectangular cross-section. The solution of the problem may be of interest 
in the problem of a shaft, press-fitted with a collar or a wheel,+ since if the 
diameter of the shaft is very large compared with the dimensions of the 
collar or the wheel, we can consider the problem two-dimensionally as one 
of a semi-infinite body compressed by an elastic plane. 

The problem of a semi-infinite body compressed by a rigid plane has 
been considered by several writers (1, 2). But the case when it is com 
pressed by an elastic plane has not yet been treated. We shall derive an 
approximate solution of the problem and examine the pressure distribu- 
tion on the contact surface of two bodies. 

2. Let us consider a semi-infinite body compressed by a rectangular 
prism, with the width 2 and the height 2h, which is compressed by 2 
vertical force P, distributed uniformly over the upper surface of the prism 
as shown in Fig. 1. The origin of the coordinates is chosen at the centre 
of the contact surface of the prism and the directions of the axes of « and 
y ave taken parallel and perpendicular to the boundary of the semi-infinite 
body, respectively. To simplify the calculation we shall not consider the con- 
dition to be satisfied on the vertical surfaces of the prism, since it can be 
assumed that the traction o, on the vertical surfaces will not appreciably 

+ The solution for the problem of a shaft press-fitted with a collar will be given 
later paper. 


[Quart. Journ. Mech. and Applied Math., Vol, IV, Pt. 3 (1951)] 


aff 
col 


whi 
ove 
diti 





shea 
surf; 
be ni 


W 


then 


mate 




















TWO-DIMENSIONAL PROBLEM OF AN ELASTIC FOUNDATION 261 


affect the distribution of o, on the contact surface.+ Then the boundary 


conditions for the prism areé 
P 
S P. cosnzx 
2 = +} stye= 
r 0 
p ‘ (1) 
: at y 2h 


- 0 
where the P, are unknown constants. In (1) we assume that 7,,, vanishes 
over the contact surface, but the assumption is far from the actual con- 


lition. since the friction between the two bodies in contact gives the 
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Is ( 
rive shearing stress. But in so far as the distribution of pressure on the contact 
istril surface is concerned, the influence of the friction between the bodies can 

ve neglected, as the writer has pointed out previously (3). 

angular | We shall assume / to be large so that sinh nzh = cosh nzh and 
dD | sinh nzh cosh nazh+-nzh = cosh?nzh: 

e pl - 4 ° . . . 
ak | then the stresses which satisfy the boundary condition (1) are approxi 
e cel nately civel hy 
ol 
Cage > 2 N7TY)COS NTL 
1-1 — . 
the ec P L 
it, can De > r 4 v(] n7y)Ccos nTx » (2) 
yrecial 

7 S naye"7" sin n7x 
rive ray 


the Appendix. 
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The displacements in question are obtained from the following relations: 
ou l ov l Cv. ou 2(1+-y 
—(o,,—V; Cy), — = —-(o,—v,¢;), cco il ( ut 
ea E, cy =«C@#d, Cx Oy E 


, 


where E, and v, denote Young’s modulus and Poisson's ratio of the prism, 
respectively. 

The vertical displacement on the contact surface can then be easily 
shown to be 


=. P 
v ~ > —cosnmx, aty Q), (4) 
7H, <1 1 


n=1 


which is independent of the distribution of load on the upper surface of 


the prism. 


At the contact surface (1 > x 1), we may expand o, and v as power 
series in 2, . : 
Co Va. x" v Ss | ed (5 
y — U n 
n=0 n=0 


en Yc. cosmn, 
hn nm 

m 0 

in the interval 1 r 1, where 
| 4(—1)" 
0 ) Om 0 Cim 9 9 
2n+1 m=" 

4n( Ly’ = =2n(2n—1) 
ae —s ——— Co_im (m By iy Bysne) 
mtr“ m*“tr- : 


Then (5) can be transformed into 
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e.. Be COS INTL 5 ac 
Uy het ino" nm 
m 0 n 0 
(h 
x f 
v YS cos marx Nbc 
y= hot n-nm 
m=0 n=0 


Inserting these relations into (1) and (4), and equating the coefficients 0 


P., we have the following relations between a, and b,: 
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The boundary conditions for the semi-infinite body are 


- 0) 
(x l and x 
0 
S ad ee 
0 
(1 
( S h pe 
n=0 


l) 
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at y = 0. (3) 


If we assume that a point on the y-axis at a distance / from the origin 


does not move vertically, the vertical displacement on the contact surface 


is expressed as (4) 


1 1—a 
2 fr l l+», [ 

. a, (7 r)log —dr ~ a, (+2) dr 
7H, } : } 7H, J 
~ 9 = @ 

/ T 
2? lI+y f 
r (r—ax)log —dr —= o,(r—ax)dr, (9) 
7H. ‘ aa rE, 4 ‘ 
0 


where E, 


infinite body 


and Vs, 
respectively 


Equ iting the values of 


. : ° 
( YS ba = \ a () 
my 7H, — } 

] 
’ 
K \ a } v)-" log? 
: } 0 
Putting 0 in (10). we have 
4 4 a 
{ 
TwH, Lu (27 1) 


trom which ( is obtained 


denote Young’s modulus and Poisson’s ratio of the semi 


in (8) and (9), we have 


aye" I< or dr 


d) 


4log/ a Vo) = Ay (10) 
rE, kat DN 
n=0 
flog] ah typ) + a, (11) 
7H, aw 2n+-] 


Differentiating both sides of (10) twice with respect to x, and putting 
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Inserting the above expressions for the b,,’s into (7), we find the following 
simultaneous equations for the a,,’s. 


~. fa lo a; 

S a > oes: Coon 0 (m Dy Diexeds (13) 
hat \Km on Ly 23—2n+ 1 

n . 1=0 * 


where k = E,/E,. 
From (13) and the first of equations (7) we can determine the a,,’s, and 
P. is obtained from the relation 


r= > a, (m ee oe 


uw -nm 


Then the stresses in the semi-infinite body are obtained from the following 
equations (5) 


* (ay— 1)e-° sin « COS at? 


oO i da 
‘ X 
0 
% - a(ay— 1)e-® sin « COS a2 
sep | S32 do 
— n° — a* 
n i = 
1 
” (ay+ 1)e-* sin « COs aa’ 
CO, P . da 
J ¥ 
" ». (14) 
. » f a(ay+ lje-® sin a Cos ax 
2 S P. ( 1)" : 3 5 da 
— A Nn? 71*— 0" 
; 1 0 
x 
T,, =P | ye-“sin asin ar da 
0 
— v7ye~% sin x Sin aa 
12S P(—1)" DA oi da 
at P aw a _ 
” 


0 
In the particular case when £, is very small compared with F,, we can 
put k = 0, and from (7) and (13), we have 
P 


1 2 aa ’ 0 ° 


When £, is very large compared with F, we also have 


b, b, ss 0, 
“ a = a 
or S ! 0, + " —_ 
Lar Dn — 1 Lat 2N— 3 
n=0 n 0 
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265 
ins we can readily show that a,, satisfies the above equations and a, is deter- 


mined from the following relation: 


r - P 
es ” . 
= Ps b< re Ao Oo! Ao 


és 0 ~ a 


Then the normal stress on the contact surface is 





ind - 
P = 1.3.5...(2n—1 P 
7 a ye “4 . ( ) yon 
— on 2 
) 77 Wb. 7 l _ 
=, v 
which agrees with the result obtained by M. Sadowsky. 

’ In the case when the two bodies are of the same material, / is unity, and 
when one of them is of steel with Young’s modulus 21,500 kg./mm.? and 
the other is of duralumin with Young’s modulus 7,220 kg./mm.?, then k is 
2-978 or 0°3358. In these cases, values of a, calculated from (7) and (13) 
are tabulated in Table 1 (unit P). 

TABLE 1 
I 
$17 77 
“rs I2Q 
+ = 154 
( +409 
l4 ar" “e 
; I 134 
79° 115 
916 1°405 
452 1°662 


] 


ind the values of b, and C caleulated from (7), (11), and (12) are tabulated 
in Table y unit P E. Vy () 3. / l. 5). 


) 


TABLE 2 
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The pressure distributions and the vertical displacements at the contact 
surface for several values of k are shown in Figs. 2 and 3. Fig. 2 shows 


by | 
[Sb 
[t- 
13+ 




















a a 


Fic. 2 





r | 


that the contact pressure between two elastic planes has a tendency to 





concentrate at both edges of the prism, and the tendency is more pro- 
nounced as the value of k increases. The values of the stress at the edges 
for several values of k are tabulated in the following table: 


> 


TABLE 3 
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act By means of the previous calculations, we readily find that the values 
WS of (o,,) 


contact pressure o, can be represented approximately by the sum of the 


Of 02 03 OF OS 06 OF 08 09 10 
ar , 95358 0 
/ 


Q2} K=00 aad / / , 


tabulated in Table 3, are valid under the assumption that the 


max? 

















Fic. 3. 


first few terms in the series of 2?”, in the closed interval 1 x l. 
Otherwise the contact pressure obtained in this way differs from the 
substantial one when x approaches the values +-1. But even in this case, 


the solution is valid at every point in the semi-infinite body except at 


the points in the close neighbourhood of the edges of the prism. 


| APPENDIX 
In our previous discussions the influence of the traction o, acting on the vertical 
| surfaces of the prism upon the displacement v on the contact surface has been neg- 
y lected. We shall find a more accurate solution, by using a numerical method, in 
A which the effect of the traction acting on the vertical surfaces is taken into con- 
sideration. Now, consider a rectangular prism of the height 2, compressed by 
| distributed loads acting on its upper and lower surfaces, as shown in Fig. 4. If we 
\ represent 1 listribution of t load by 


Ss P. cosn7er, 
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The stresses which satisfy the boundary condition (15) are (6) 
o pa P’ & (y)cosn7x 4 Q, Ny(X)cos n7ry} 
n=1 | 
| . 
P , : - (16) 
o = > {Pi nny)cos nx + Q, &,(x)cos ny} | 
ka 1 | 
whe! 
sinh nz cosh nay + ny sinh nz sinh nay — nz cosh nz cosh nay 
j : _ : 
sinh nz cosh nz + nz 
sinh nz cosh nay — nzy sinh nz sinh nay + nz cosh nz cosh nzry 
j . — > 
sinh nz cosh nz +-n7 
and P’ are the solution f the simultaneous equations 
I > 2 > 1)2I(7,n)L(m, 7) (n yy wet (17) 
l 1 
/ 4(— 1)/nj? sinh®nz 
in whicl . ————s . 
7(n* +-9*)*(sinh nz cosh na + nz) 

Using (16) in place of (2) and the approximations sinhnz = coshnz and sinh nz 
cosh na --na7 cosh*n7, the more accurate formula for v on the contact surface is 
the con nt term being omitted) 

, 
5) ~ P’ 
= S = cosnaz, (18) 
TE, ya) n 
n 1 
since the displacements v due to the stress a, and the second term in the expression 
Io val h on the contact surface. 
Accordingly, in our accurate calculation, (13) should be replaced by 
" 1 a l ne a 
* S j a 0, (19) 
9; 9 nm 
homme hint VO Lam = oll l 
r l }=0 
where a Po ir. 
\s an example, we shall consider the case for k 1, and take the values of a, 


obtained in the previous example and the corresponding values of P,, as the first 
approximation. Inserting the numerical values of P, into (17), and solving the 


simultar is equations, we can find P? 


and subsequently the corresponding values 


of x, for the first approximation, as are shown in Table 4. Inserting these values 


f ~,, int 19), and solving the simultaneous equations, we have the values of a, 
for the second approximation. Repeating this operation, we have finally the values 
of x, which do not vary by the operation, and subsequently the corresponding values 


of a,, which are the ones for the exact solution. 

The p1 ire distribution on the contact surface, calculated in this way, is shown 
in Fig. 5+ and the approximate one, previously obtained, is also shown by a dotted 
line in the figure. From the figure we see that the influence of the normal traction 

ting on the vertical surfaces of the prism, upon the pressure distribution on the 


contact surface, is not perceptible except at the edges of the prism. The accurate 


The pressure distribution in the second approximation almost agrees with that for the 
nation. 
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TABLE 4 









































Ist approximation ap 0°417, (Oy) max = —1°26P 
, I 2 3 4 5 6 | 7 8 
a | ot ~f - 
an orr52 0°634 sO} 2°006 O°ol7 | 2°062 | O79 | ogI6 | 0°452 
Nos 1°25 I*4I } 1°49 1°55 | 1°57 1°58 1'60 | 1°64 
2nd approximation @y = —0°423, (0,)max riiP 
hi I | . 3 4 5 6 | 7 8 
ene ——_—_—_| : a - -| st 
An or125 0052 0556 2°d71 2°360 0°360 | 1°d6¢6 2°353 
Ney 1°24 | 1°43 I°52 1°60 1°64 | 1°65 1°69 1'69 
3rd approximation dp 0497, 0). 1:10P 
| | | im ; 
, I 2 2 | 4 | 5 6 ' | 8 
= a — See = | _—— 
| . 
an OrrIs o451 2°507 4°975 | 27051 | o080 4°330 4°168 
hes 1°24 1°43 1°52 | 1°60 1°67 | 1:69 | 1°70 1°70 
value of (Gy)max 18 —1:10P, while the approximate value previously obtained is 


1-26P, which is about 15 per cent. larger compared with that obtained here. 

From these results we can conclude that the approximate solution previously 
obtained, where the effect of ¢, acting on the vertical surfaces of the prism is neglected, 
can be used with sufficient accuracy for the calculations of technical problems. 

In conclusion the writer wishes to express his thanks to Miss E. Itagaki, his 
assistant, for her earnest help in numerical calculations. 
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THE MOTION OF AN ELASTICO-VISCOUS 
LIQUID CONTAINED BETWEEN COAXIAL 
CYLINDERS. | 


By J. G. OLDROYD 
Courtaulds Limited, Research Laboratory, Maidenhead, Berks.) 
Received 20 December 1950] 


SUMMARY 
The theory of periodic and steady states of motion of a liquid in a narrow annular 
gap between long vertical coaxial cylinders is considered, when the viscous and 
elastic properties of the material at small rates of shear are measured by three con- 
stants, a viscosity coefficient and two relaxation times. When the outer cylinder 
undergoes forced harmonic angular oscillations about its axis, and the inner cylinder 
is constrained by a torsion wire, the amplitude of uniform oscillation attained by the 





inner cylinder is calculated. The ratio of the angular amplitudes of the two cylinders 
has been plotted as a function of the frequency, for different values of the constants 


of the material and of an available apparatus. From observations of the frequency, 


amplitude-ratio curves for some actual liquids which are approximately of the type 
onsidered, values can be deduced for their elastic and viscous constants. In steady 


flow between rotating cylinders, elastico-viscous liquids behave like purely viscous 


iquids whose viscosity is a function of the rate of shear. Additional information 
about the rheological equations of state of actual liquids can be deduced from steady 
rotational experiments. 


1. Introduction 
Iv is well known that elastic properties of liquids in shear can be detected 
and measured by observing suitable types of oscillatory motion (see, for 
example, 1). For the purpose of investigating the rheological properties 
of dilute polymer solutions (2), an elastoviscometer has been built in this 
laboratory which provides a means of maintaining a liquid in a state of 
uniform oscillation in a narrow annular gap between vertical coaxial 
cylinders (3). The outer cylinder can be either rotated steadily at any 
speed between 0-3 and 3,000 r.p.m., or maintained in a state of harmonic 
angular oscillation about its axis at any frequency between 0-25 and 25 c/s; 
either the constant deflexion or the amplitude of uniform oscillation at- 
| tained by the inner cylinder, when constrained by a torsion wire of known 
restoring constant, can be observed. It is the purpose of the present paper 


to outline the mathematical theory underlying the quantitative interpreta- 





tion of observations made on different liquids. Only two-dimensional 
motion need be considered, since it is known that end-effects in steady 


(Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 3 (1951)] 
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rotational experiments can be eliminated by repeating the experiments 





in all cases of physical interest. 









































f 


with the apparatus filled to different heights; in oscillatory experiments, al 
in which the annular gap is always completely filled, the observations on ae 
Newtonian liquids of known viscosity have confirmed that end-effects are pro 
negligible. 
= . : . in § 
It is hoped to establish the rheological properties of certain polymer pro 
solutions in two or more stages, the first stage being the investigation of pe 
the behaviour in a simple shearing motion at small rates of strain. It has hen 
been found that this behaviour can be approximately described by a linear oti 
equation > LA, > nly AY) (1) prec 
invé 
relating the velocity gradient y, the associated shear stress 7, and their | y 
rates of change y and. The properties at sufficiently small rates of strain oft 
are measured by three constants, a coefficient of viscosity 9, and two po 
relaxation times A, and A,, with the direct physical significance that if the -_ 
motion is suddenly stopped the shear stress decays as exp(—t/A,), and if | und 
the stress is removed the rate of strain decays as exp(—t/A,). The viscosity | diff 
No is that obtained from experiments in the elastoviscometer involving | on 
slow steady rotation of the outer cylinder. The values of A, and A, can be —_ 
deduced from observations of the amplitude of uniform oscillation of the 4 ..., 
inner cylinder, plotted as a function of frequency in the range 0-25-25 c's, diet 
when the outer cylinder undergoes forced harmonic oscillations of known 
small amplitude. The estimation of A, and A, is made by comparing an 2.1 
observed frequency-amplitude curve for a given material with a set of W 
theoretical curves based on the equation (1) with different values of A, and 4 the: 
A,. The justification for associating the equation (1) with an actual material inde 
depends on the results of experiments with different inner cylinders and diste 
torsion wires being consistent with the same values of 7p, A,, and Aj. tion 
An equation of the type (1) was first proposed as a possible stress, rate- 
of-strain relation by Jeffreys (4). It was shown to be associated with a It is 
suspension of elastic particles in a Newtonian liquid by Frohlich and Sack nd 
(5), and some of the consequences of such an equation have been discussed 
by Burgers (6). With (1) must be associated the auxiliary condition that 
A, 7—7pAsy is continuous at a discontinuous change in stress (cf. 5). It ina 
follows that, at a sudden application of a shear stress to the material at bebe 
rest, the initial rate of shear is A,/A, times the corresponding rate in a steady 
state under the same conditions. This factor can only represent an increased 
response to suddenly applied stress due to elasticity, so that whic 
Ay : As (2) mae 
would 
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A single equation (1), associated with simple shearing of a material at 
sufficiently small rates of strain, together with the usual assumption of 
negligible bulk compressibility, does not suffice to specify the rheological 
properties of the material completely. In order to describe the behaviour 
in general motion, a set of equations of state with the right invariance 
properties is required (7), and there is not a unique way of recasting (1) in 
an invariant form. On the other hand, in order to present a mathematical 
treatment of the type of motion involved in the determination of A, and A,, 
which is not strictly a simple shearing motion, it is desirable to specify 
precisely the idealized material to be studied and therefore to select an 
invariant generalization of (1) on which to base the theory. 

When it has been established that an equation (1), with certain values 
of the constants, can be associated with a particular liquid, one must then 
consider what further experiments are required in order to distinguish 
between the various possible sets of equations of state which reduce to (1) 
under simplifying conditions. In the present paper it will be shown how 
lifferent generalizations of (1) can give rise to different types of non- 
Newtonian viscous behaviour in steady-state experiments.t Hence the 
investigation of rheological properties might in principle be carried to a 
second stage by making quantitative observations on the liquid in the 


elastoviscometer when the outer cylinder is rotated steadily. 


2. Idealized types of elastico-viscous liquid 
We consider (almost) incompressible elastico-viscous liquids in which 
the stresses p,, at any point may be considered as the superposition of two 
ndependent stress systems, one represented by a tensor p;;, related to the 
listortion, and the other an isotropic tensor —p’g,, related to the dilata- 
tion F - P 
Pik—P Gir: (3) 
lt is of interest here to examine those in which the relation between pi, 
ind the rate-of-strain tensor e\) reduces to 
ith, — he 2no( I Ay Sele (4) 
“a I’ - Joe if” k 
+] 0 E t, i 
a sufficiently slow simple shearing motion, and in which the relationship 


tween the scalar p” and the dilatation reduces effectively to 


O (all p”), (5) 

vhich represents ideal incompressibility. The notation is that of the earlier 

The possibility of qualitative differences in behaviour has already been demonstrated (7) 

stablishing that some generalizations would show the ‘Weissenberg effect’ while others 

ud not. We are concerned here with differences which can be measured by recording the 
iple on the inner cylinder at different speeds of the outer. 


T 
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paper (7), in which liquids of this type were used for the purposes of illustra- 
tion: g,;, denotes the metric tensor of a general curvilinear coordinate 
system x fixed in space, covariant tensor suffixes are written below, contra- 
variant suffixes above, and the usual summation convention is assumed. 

In order to convert equations (3)—(5) into a possible set of equations of 
state of general validity, it is only necessary to recast (4) in an invariant 
form which reduces approximately to (4)—and therefore implies (1)—in a 
sufficiently slow simple shearing motion. This can be done in an infinite 
number of ways (cf. 7). When the incompressibility condition (5) is taken 
into account explicitly, the distinct invariant generalizations which are 
linear in the stresses, and quadratic in the stresses and the space deriva- 


tives of velocity taken together, are 


D\_, dD 
1+A Din. 21 (1 LA, —]ep, (6)+ 
1 a k lo 2 > ik )) 


and equations obtained from this by adding an arbitrary linear combina- 
1 > AL) p(1) (1) 
tion of Cim ©k wit im 


which each side is a symmetric tensor are relevant, since the required 


py”, and e) p?”" to one side. Of these only equations in 


number of independent equations relating stress and rate-of-strain compo 
nents is six. Hence the simplest} generalized forms of (4) are 

(1-41 2) pie 2eu(AR rit+X, pm) = 2na(1+As SA —Snoaeitem 

(7) 

where «,, K, are arbitrary scalar constants. More complicated equations, 
involving contracted triple and higher products of the stress and rate-of- 
strain tensors (and, in general, arbitrary functions of scalar products), give 
other distinct generalizations. Such equations might well have to be in- 
voked in order to explain the behaviour of actual materials, but the set of 
equations of state (3), (5), (7) suffices to illustrate how different types ol 
behaviour can arise in steady-state experiments from different generaliza- 
tions of (1). 

The particular values x, = A,, kK. = Ay define the idealized material B 
of the previous paper (7), which was shown to exhibit the ‘Weissenberg 
effect’ in steady rotation; they permit (7) to be written, after raising 


D | D 
1+ A, ; p'tk 2y{ 1-4 As : Dik, (3) 
dt “Dl 


The particular values «, 0, Ky 0 define the material A which would 


suffixes, as 


+ D/dt denotes convected differentiation of a tensor quantity in relation to the material it 
motion, as defined in (7). 
There is no invariant generalization which is linear in the space derivatives of velocity 
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not show this effect. These two materials are distinguished from the other 
possible types because, in either A or B, a simple shearing motion is 
possible in which (1) is satisfied exactly, without restriction to small rates 
of shear. More precisely, a velocity distribution given by v; = [0, v(x, t), 0], 
referred to orthogonal Cartesian coordinates x', can be maintained by non- 
vanishing stresses pj,(21,?), p4,(21,¢t) alone in A—or p},(a1,t), p5o(a?,t) in 
B—where : i 
Py \, = lo ~ be oe 
~ ct oa} ~ ctea! 
without approximation. It is seen that either of these types includes, as 
a limiting special case when A, and A, tend to equality, the incompressible 
Newtonian liquid in which 
» l h » 


Pi 270 (and, equally, p’ 2c). 


3. Flow between coaxial cylinders 
A liquid defined by the equations of state (3), (5), (7) is now supposed to be 
contained between long vertical coaxial cylinders. All tensor quantities 


associated with the liquid are most conveniently represented by their 


physical components in cylindrical polar coordinate directions r, ¢, 2; z is 
measured vertically upwards and r = r,, r = ry (> r,) represent the rigid 
boundaries. ‘Two-dimensional axially symmetric flow is assumed, with 


velocity components 0, rw(7,t), 0. The relations between w and the physical 


components rd’, 17”, dd’, 22’ of pi; all functions of r and ¢ only, are obtained 
in the same way as for steady flow of particular types of liquid (cf. 7, 
pp. 539-40), except that partial time derivatives cannot here be omitted. 


These relations are, from (7 


f —~ 2 
CTO CW , —~, CW Cc°W 
} Ay— kK, lbh K,7Tr | no! ( LX, — H (9) 
ct cr cr ” € tc a 
crr Cai}; CW : 
} A, 2(A, Co rd 2 no(Ao K)(r . (10) 
l i 2 2 - 
ct cr co? 
— Chd ' Cw, “ Cw\? 
bad A, ; “Ky, rd =No Koll } - (11) 
¢ or cr 
o oar » 
A, czz /ot Q. (12) 


The equations of motion, expressed in terms of the same quantities, take 
the form 


orr ry dd cp ' crd’ = 2rd’ Cw 022’ c p’ 
prw-, pr ; t+ pg. 
cr } cr cr r ct C2 Cz 


(13) 


where p is the density of the liquid and g the acceleration due to gravity. 
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The observations which can be made in the elastoviscometer can give, 
at the most, w and rd’ at the boundaries, so that the equations of imme- 
diate interest are those not involving rr’, bd’, 22’, or p”. Inthe steady state, 
the relevant equations are the second of equations (13) and the following 
equation, obtained by eliminating 77’ and $d’ from (9)-(11) after putting 
C/ct 0: 


ae 1 +-2(Ay Ko +k, AQp— 2k, Kg)(% Cw cr)? éw 
2 2 2 , 


rd’ No (14) 


l+-4«,(A,—k,)(7 6w/er)* or 
Hence the observed behaviour, so far as is revealed by measurements of 
the torque on the inner cylinder at various speeds of the outer, will be 
precisely that of a non-Newtonian viscous liquid whose viscosity 7 is a 


simple function of the velocity gradient y defined by 


[1+-4«,(A,—«,)y? |]. (15) 





7) Mol 1+ 2(Ay a+ Ky Ay— 2K Key)y” 
The observed viscosity is the constant y, only in special cases, for example 
in a material of type A or B. For general values of «, and Ky, mo is the 
limiting viscosity at low rates of shear. When the annular gap between the 
cylinders is sufficiently small (or, when 7 is a constant, without this restric- 
tion) the usual formulae for a Couette viscometer, 
27% Q (r3—rj)G 


7 . (16) 
hQ 


4rré : 


non 


define the y, 7 relationship parametrically. Here, 2 denotes the angular 
velocity of the outer cylinder, G the couple on the inner, and / is the height 
of liquid between the cylinders plus an end-correction; this correction is 
readily eliminated by measuring the couples at different heights of filling. 
In this way different values of «, and «x, could be distinguished experi 
mentally if A,, A,, and yn») were known. 

In the case of periodic motion of frequency », it is observed that there 
is not an exact solution of equations (9)-(13) in which w, rd’, rr’, }¢’ all 
vary sinusoidally, because the equations are not linear in these quantities. 
For small rates of strain, however, there will be a solution in which 


_ 


rd’ O(nor Cw/er), rr’ Oj No A(T Cw/é r)?}, dd’ O| NoA(r Cw ¢ r)?), 


ar c rd’ CW c 2w » 
and rp’ +A, —— = r( bite —— (17) 
ct or ” oter 


9 


with an error of the order of 7)A3?(r éw/ér)®. Equation (17) is exact in the 


case of a material A or B. since in A dd’ is seen to vanish, and in B rr 


vanishes. The following theory is an exact treatment for either of these 





prototypes, and will hold approximately when «, and «x, have other values. 
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Equations (13) and (17) permit a solution of the form 
rd’ re f(rje27rt, w re F(r)e27'nt, (18) 
where f(r) and F(r) incorporate appropriate phase factors. The equations 


relating rd’ with w reduce to 


d(r*f )/dr 2rinpre F (19) 
and (1-+-2m2nd,)f nol 14-27ind,)r dF dr. (20) 
Eliminating F’, we have 
d? d : 
re 7 | q2r2— 4] f 0, (21) 
at ar } 


where a(n) is defined by 


» “4 » “ 
: 27inp( 1 2r7ina ) 
x" . : . ( har < arga < }n). (22) 
Holt 27indAs) 
The solution of (19) and (20) can therefore be expressed, in terms of Bessel 


functions of complex argument, as 
f aJ,(ar)+-bY,(ar), F x| ad, ( wr) + bY; ( vw’) | (27inpr), 23) 


where a and } are complex constants to be determined by the boundary 
conditions 

If 6,, 8, are the amplitudes of angular oscillation of the inner and outer 
cylinders, J is the moment of inertia of the inner cylinder about its axis, 
L the length of the annular gap, and K the restoring constant of the torsion 


wire, we have 


F(r,) 27ind,, F(r.) 27nd, e*, (24) 
where ¢ is a real constant. Also, from the equation of motion of the inner 
cylinder, 2ar, Lf (r,) (K —4n?n71)6,. (25) 


These equations determine a, 6, c, and 6,/6, in terms of the frequency. The 
most readily measurable quantity in oscillatory experiments with the 


elastoviscometer is the amplitude ratio 





27 ) 7 Vs L - 7 
H salt Sl | Ja( xr", )¥4 (a7) Y5(cxr'; Jy (axr’,) | 
Un) l - 7] x | 26) 
| —! -| J, X JY (ar's) Y,(a7,)J, (ar) | | 
[ K . 
| a l [ J (ar, )Y4 (ary) Y,(a7,)J, (a7) | 


It is observed that a radial variation of vertical normal stress zz 


p”) is required to maintain strictly two-dimensional motion. It is 
easily verified that, in addition to the variation arising from the centripetal 
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acceleration prw?, the variation of zz across the annulus is of the order of 
oA, (%2—1,)r(@w/éer)*. The departure from two-dimensional motion which 
must arise when a free surface is present at the top of the annular gap can 
therefore be made negligible by sufficiently reducing the annular gap and 
the amplitude of oscillation. In practice, a check is made experimentally 
that the amplitude is sufficiently small for the above theory to hold, by 
noting that the amplitude ratio is (very nearly) independent of the absolute 
value of the amplitude of oscillation of the outer cylinder. 


4. Numerical calculations 
The numerical constants associated with the three different inner cylin- 


ders at present available for oscillatory experiments are as follows. 


L (em.) | I (gm. cm.*) 


Cvlir er | ry (cm.) Vo (C71.) 
] 2°1212 2°1712 25°40 | 2677°5 
II 2°0701 2°1712 25°40 =} 2468 
19692 2°1712 25°4 


"40 2095 





Calculations have been made with two values of K, corresponding to the 
torsion wires which have been used in most of the experimental work of 
Toms and Strawbridge (2), namely 0-171 x 10° dyne em./radian (Wire 5a) 
and 2-79 x 106 dyne cm./radian (Wire 8). 

The exact formula (26) is unsuitable for direct computation because there 
appear to be no tables of the Bessel functions J, and Y, for general complex 
arguments. An approximate formula, valid for sufficiently large n, is ob- 
tained by substituting for the Bessel functions their asymptotic expansions 
in powers of (ar)-!. After some reduction, the following approximate 


formula is obtained 


l 3(r.—r.) 8) 
(1 ",) Joos(r ro 


8r2 
15r.—3r 15 3 45 y , 
_ieek ee satl(*tsn-- ass S|sin(r,—7,)a|, (27) 
Sre ax L678 ax? Sr2a 128r§a*] _ . 


where 7) = 4/(7,72) and 


S = L K ae . (28) 
2apr3 L\422n? 


When the above numerical constants are used, it is estimated that the 
formula (27) gives # with an error of less than about 1 per cent. for values 


of n sufficiently large to make |a| > 2 em,~! 
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It follows from (26) that the slope of the (n,#) curve at the origin is 


given by 


. O(n) Sarre rs Lr 
lim be St 10 (29) 
no 7 (r3—rj)A 


a result which is independent of A, and A,. The observation of this slope 
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Fic. 1. The variation of #(n) with A, and As when p l. No ° 8 


(Inner Cylinder IT, Torsion Wire 8). 


in an experimental (n, #) curve gives an independent, though less accurate, 
method of obtaining ». (It has been found that the values of 7, so ob- 
tained for polymer solutions are consistent with those obtained from 
steady-state experiments 

The function 3(n) depends on six parameters which can take different 
values when the liquid or elastoviscometer components are changed: p, 7p, 
Ay, Ay, , and A. An example of the variation of the theoretical (n, #) curve 
with 7, and A for a given liquid has already been quoted by Oldroyd, 
Strawbridge, and Toms (3, Fig. 7), and this need not be considered 
further. The results of numerical calculation show that the form of the 
curve is very little dependent on p, when the other parameters are fixed, 
in the range 0°75 p 1-5 (g./ml.). The manner in which the curve varies 
with A, and A, is illustrated in Fig. 1, and the effect of changing ny, other 
things being equal, is illustrated in Figs. 2 and 3. (Numerical values 
marked on the diagrams are in C.G.S. units.). The curves for Newtonian 
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liquids (A, = 0, A, = 0) show a peak in #, not exceeding unity, at a value 
of n almost exactly equal to the frequency n, of free oscillation of the inner 





























cylinder and torsion wire (defined by S = 0). Those based on values of 
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Fic. 2. The variation of #(n) with ym» in the case of Newtonian 
liquids (A, 0, A, 0, Inner Cylinder II, Torsion Wire 8). 
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Fic. 3. The variation of #(n) with ny when A, = 0-05, 
A, = 0-01 (Inner Cylinder IT, Torsion Wire 8). 
r, 0 show a higher peak at a frequency greater than n». The height o! 


the peak increases with A, and decreases with increasing A,. 

The change in shape of the (n, #) curve with A, and A,, at constant density 
and viscosity (cf. Fig. 1) is sufficiently marked to permit approximat 
values of A, and A, to be assigned without ambiguity to an experimental 


curve, by comparing it with the theoretical curves based on the known 
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and 7», and a series of values of A, and A,. For example, the values 
\, = 0-065 sec., A, = 0-015 sec. were assigned in this way to the solution 
of poly(methy! methacrylate) in pyridine (p = 0-98 g./ml., 79 = 7-9 poises) 
the curves for which were given by Oldroyd, Strawbridge, and Toms (3, 
Fig. 6); the observed curves for the different cylinders and torsion wires 
were in general agreement with the predicted curves based on these values. 
A device which has been found useful in assigning values of the relaxation 
times rapidly is a graphical method such as the following. Two frequencies 
n, and n, (> mp) are chosen, at which the variation of # with A, and A, is 
considerable, and contours A, = constant, A, = constant are constructed, 
on a plot of 3(n,) against #(n,), from the theoretical curves. Approximate 
values of A, and A, can then be assigned immediately by marking the point 
H(n,), Hn) on this plot which corresponds to an actual material. The 
choice of n, and n, depends on the viscosity and density of the solution 
being examined. 

In order to avoid extensive calculations for the analysis of every new 
experimental curve it has been decided to study polymer solutions of such 


compositions as to give 2” poises only (m 0, +1, +2,...). The 


systematic computation of (n,#) curves has been carried out for these 
viscosities and for p | g. ml. (A small correction is necessary before 
these are suitable for direct comparison with experimental results for 
liquids of appreciably different density.) It has been found sufficient to 
consider values of A,, A, increasing from zero in steps of 0-01 sec. up to 
about 0-08 sec., in order to assign values to the experimental curves so far 
considered. Most calculations have been done to correspond to the experi- 
mental set-up of Inner Cylinder II and Torsion Wire 8, and the first com- 
parisons of observed curves have been made with the resulting theoretical 
curves. 

The experimental results from steady rotational experiments show that 
the liquids exhibit marked non-Newtonian behaviour (ef. 3, Fig. 5). The 
variation of » with y, calculated from equations (16), indicates that more 
complicated equations of state than (7) are required to explain the rheo- 


logical properties of the liquids so far observed at finite rates of shear. 
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THE SURFACE ELEVATION IN CELLULAR 
CONVECTION 


By HAROLD JEFFREYS (St. John’s College, Cambridge) 
Received 6 June 1950] 
SUMMARY 

3énard observed that in cellular convection in a fluid heated below, the surface 

was depressed over the rising currents. His theoretical explanation is fallacious, and 
an analysis leads to the opposite result. A discrepancy remains between theory and 
observation. An alternative method of numerical solution was used as a check and 
may be capable of extension to other problems. 
1. Tue formation of cellular convection currents in a fluid heated below 
was studied experimentally by H. Bénard (1) and there have been many 
theoretical discussions. My attention was recently called by Miss M. E. M. 
Bland to a passage (1, p. 92) where Bénard states that the currents rise 
in the centres of the cells and that the surface of each cell is concave 
upward; thus to each centre of ascent corresponds a depression of the free 
surface. This appeared surprising in the light of several problems that 
have been solved with a view to the explanation of thermally produced 
winds. In these the disturbance of temperature is taken to be a function 
of the horizontal coordinates, varying slowly with the time, and the result 
in each case is a slow motion with the pressure in a heated column increased 
at the top and diminished on the bottom; while the flow at the top is out- 
ward from such a column. Thus the flow at the top is away from the 
elevated parts. 

[t appeared that the different conditions of the problems might provide 
an explanation, but a thermodynamical argument suggests that they 
cannot. In the steady state the dissipation of energy must be compensated 
by a supply, which can come only from the contraction of elements of the 
fluid where the pressure is high and expansion where it is low; that is, for 
an element describing a circuit, { p dv must be negative. Since the fluid 
is assumed incompressible, and since the temperature is taken constant 
over the free surface and the bottom, this implies that a particle must be 
under excess pressure when it is moving towards a place of lower tempera- 
ture, that is, upwards; and hence the elevations must be over the rising 
currents. 

Bénard gives a qualitative argument to show that his result could be 
inferred on physical grounds, by considering the currents near the bottom. 
These are being driven against viscosity and imply a decrease of pressure 
along the bottom from a region of descent to one of ascent. But the 


{Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 3 (1951)] 
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conclusion does not follow for several reasons. (1) In the corresponding 
problems of wind formation the pressure gradients along level surfaces 
have opposite signs at the top and bottom, and this might be true here 
too. (2) The relevant stress components for bottom currents would be 
P11» Po, in a usual notation; while the surface elevation determines p,,, 
and the differences, even at the same place, are not negligible in a motion 
largely controlled by viscosity. (3) The shear stresses p,3, Po, are not zero, 
and play a part in driving the bottom current. Unfortunately (2) also 
vitiates the thermodynamical argument, because the normal stress at the 
free surface of a viscous fluid would differ from —p. 

It might be argued directly that the rising current can be stopped at 
the free surface only by an opposing vertical stress, which can be provided 
only by an elevation of the surface. A full discussion is desirable. 

2. Notation 
Xo, Xz = xX, y, z coordinates; z vertically upward. 

Uz, Ur, Uz = U, Vv, w Velocity components. 

Pi = Stress components. 

p = density. 

x = thermometric conductivity. 

v = kinematic viscosity. 

|)” = temperature. 

8 = eV /éz in undisturbed state. 

x = coefficient of volume expansion. 

h = depth of fluid. 

Accents refer to departures from the undisturbed state. 

Departures are taken to be of the forms (cos/x, cos my)(sin la, sin my). 

It was only in the work of Pellew and Southwell (2) that the solution 
was first found for the motion with a free surface, and no discussion has 
attended explicitly to the variations of pressure and of surface elevation. 
I use the method of my original paper (3) to state the equations; the formal 
solution is found by Goldstein’s method (4, section 14.062). This allows 
for Pellew and Southwell’s adaptation to a free surface. It is necessary to 
notice, however, that in all treatments of the problem the equations ol 
motion have been taken in the form 
: (u,v, uw) (- ‘ ; , ; )p + pvV2(u, v, w)— (0, 0, gp), (1 
Ox Cy CZ 
which implies a special and rather unusual definition of pressure in relation 


to the stress components; (1) in fact is true only if 
ae (Ou, . Oru, Cu x 
Pik —( p-- pv@)o;;. mal .t. - ‘) G = =. (3) 


CX}. Cc v; 
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\ding Since the fluid is taken to be incompressible, this is permissible because 
faces the pressure does not appear in the equation of state and is required only 
hones to be a scalar. But in studying the surface elevation, the only datum will 
ld be be the variations of normal stress over z h, and we need the relation 
3S —p between p and ps3. Considering only disturbances from the static state we 
siti have, with ¢@/ et 0, 
Zer Su KV=)"", 6 KV2AV"” (3) 
ils 
Cp K OV 7 sa 
t the ~ r V4i gapt - i). (4) 
Un 4 
c= = = oV’ i 
ed a Vp ve pvV 4 Jxp - 0, (9) 
vided oz 
‘ KpV ¢ * id : 
whence Va p : = V44 - vk pvV 4] 0. (6) 
With [2 +-m*)h? a*, (7) 
Pimn 2 
ve have p ~ Vip’: (3S) 
Si 
und from (2), (3), (6) 
[ KpV ¢ - ry 
Dp fe Ved | vKpvV *] 
i” 5 C2 
2pvK ¢ — 
KV 2] : V2) (9) 
») c 
p ¢ ([h> ei vP 
—- V4f 2V2V'l, (10) 
> oz \a* 
my . . ° ° . 
| ti since the second and third terms of (9) can be neglected in comparison with 
uth . . ; 
the first and fourth when ah is small. Put now 
n | 
Op 4 
sai ; gaph 
ati a th 2 hé 7. pL - I r (11) 
orn KI , 
illo s 1 d? i 
: Then Ps WV | —- b*\( aes abe)", (12) 
uT'y 1E \aE- dg" 
ms 4 Pols > 
where M kvp7® Bf *h® = gah/b*zp: (13) 
and therefore M, are positive since f is negative. 
Take : . ; 
| f(a, y)Z(z {(x.y) > A,, sin n€. (14) 
] + 
ath a 
Then 
Das M|—3B-+- > {(n?+-b?)(n?+-3b?)nA, — Bheos n€]. (15) 


. * In (10) of (4) p. 414, 7 should be 2. 
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Since the case of a free surface corresponds to sums over even 7, we use 


(18) of (4) and find 


| * b2(n4— b*)-+- pb? - 
Pas MB S ———_ f ~ cos n€}. (16) 
2 <4 (n?+b?)8—pb? 
Put n = 2m, b = 2s; since we are concerned chiefly with € = 37 we need 
also the sums 
“ ( ] )" ] = 
5 - = - + cosech zs. (17) 
a m* +s" 28s“ 28 
1 
~ (—l1)" I 7 a 
> - y ; +... cosech sz + —. cosech sz coth sz. (18) 
<= (m*--8*)° 2e* 88 4s° 
Then 
Pss M B{ 3x8" cosech sz coth sz 
Moy? 2) 1 102 
X (—1ym m*(m=-+- 3s") us 
— (m?-+-s)*{(m?-+- 8s?) — | us*! 
(19) 
With the numerical values (4. p. 415) 
b 1-70. pl LSO-S. (20) 
this gives 
P33, MB 0-501 4- (2-854-+ 0-0929— 0-0099-L ...) +2436 (21) 


in which a correction has been made for the probable value of the remain- 


der. Also A. 0-01024B (22) 


and sin 2€ is the only harmonic that does not change sign between 0 and }7. 
Hence a column of heated fluid corresponds to negative b, therefore to 
negative ps, at z = h, and therefore to an elevation of the free surface. 
Further, V2V’’ has the opposite sign to V’, and 8 is negative: hence w is 
positive in the heated columns, and there is an elevation of the free surface 
over the rising currents. as expected. Thus there is an unexplained dis- 
crepancy between theory and observation. The surface elevation ¢ satisfies 


y 


¢ Ja - 9: 

. P33 __ (15h (23) 
A, gp 

nearly. It must be remembered that / in this formula is twice the depth. 

Since « is the coefficient of volume expansion, the surface elevation 1s 


comparable with the direct effect of linear expansion. 


In view of the result and of the possibility of a mistake in sign in the 


calculation, I sought for another method that would serve as a check. 





The formal solution of Pellew and Southwell would also, I think, involve 
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a risk of errors in sign, and the safest method appeared to be a numerical 
one. I first tried the method of finite differences of my original paper. 
This had never been applied with the correct boundary conditions. A test 
of the applicability of this method was desirable for other reasons. Both 
the Fourier series and the exponential solution are suited only for a 
differential equation with constant coefficients, and more complex cases 
may need to be considered; but a numerical method should still be 
applicable. We now take € = z/h = w€,A = v4. The method was applied to 
the case of two conducting rigid boundaries and to the lowest (symmetrical) 
mode. With intervals $, A was found to be 604; with 4, 955; with 4, 1180. 
Extrapolation to zero interval by Richardson’s method gave A 1548. 
As the correct value is 1708, the convergence is far too slow to be satis- 
factory. A possible explanation is that to form a centred third difference 
it the boundary it was necessary to take two extrapolated values of V’; 
but in fact only first and second centred differences were used in forming 
the boundary values of the first and second derivatives. With two extra- 
polated values available it would be possible to allow for third and fourth 
differences. When this was done, however, it was found that the signs of 
the extrapolated values were reversed. Thus it appeared that the method 
is precarious when used for a differential equation of high order containing 
a large numerical coefficient 

\ method briefly mentioned in my first paper, but not completed, was 
then tried. This was to make up a polynomial that satisfies the boundary 
conditions exactly, and then to choose A so as to satisfy the differential 
equation where the solution is a maximum. Since there are three con- 
ditions* Z = 0, Z 0, Z” —a®Z’ = 0 at each boundary, and since the 
solution must be symmetrical or antisymmetrical about the mid-point, it 
seemed likely that the solution is very fully determined by these conditions 


alone. For the symmetrical case the simplest polynomial that satisfies 
them is . . ' 9 +\92 2) 7: ¥\s 
Z ¢(1—f)+-27(1 —Z)?+- (2+ 4a?)f3(1—Z)?. (24) 


Then A is chosen so that the differential equation 


J2 3 
| a2) Z AazZ (25) 
is satisfied at ¢ | For various values of a the results were 
a> S 9 10 
A 1S70 1833 1888 


and by interpolation the minimum is for 
a 8-90, a 2-98, A 1833. 


Here accents denote derivatives. 
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As the correct values are a = 3-14, A 1708, these results are very satis- 

factory for a first approximation. It could be improved by including an 

extra term in Z and fitting the differential equation at another point. 
For the antisymmetrical case a suitable polynomial is 


Z = (1—2¢)0(1—f){14+- 32(1—Z) + 10f7(1—Z)?++- 4a2l2(1—Z)?}. (26) 
and by fitting the differential equation at ¢ = } we get 
a? 15 20 22 2% 25 
A = 17,539 16,839 16,772 16,796 16,891 
By interpolation the minimum comes at a? = 21-7,a = 4:66, A = 16,771, 
The correct values are a = 5-34, A 17,610. Again the method leads 


easily to a fairly good approximation; A must be divided by 16 to give the 


value for depth h with a free surface. 


Returning now to (12) we have to evaluate at € = 3 
Z—4a2Z" +-3a4Z’ 2520 —450a2—4-3la!—a®/64 14745 
with @* = 21-7. At ¢ 1 which is near the maximum of Z, we have 


3 
Z = 0:25. The sign is therefore confirmed. We also have, nearly, 


ss 0-16ah 
9p V max) 
in good agreement with the Fourier series solution. 

Apart from its use as a check in the present problem, this method may 
be capable of extension, particularly if generalizations of the stationary 
principle proved by Southwell and Pellew should be possible. It seems 
that such generalizations should exist because the required value of A is 
the lower bound of values such that distributions exist that make the 
supply of energy to the eddies enough to balance the loss by dissipation. 

It appears desirable that Bénard’s experimental work should be repeated. 
Chandra (5) has made a study of the conditions in air with rigid boundaries 


above and below, and finds that the motion is downward in the centres of 


the cells, whereas with a free surface above Bénard found it upward. He 
also found a second type of convection, which is not covered by the present 
theory. Similar work in CO, has been done by Dassanayake: for informa- 
tion about this I am indebted to Professor O. G. Sutton. 


REFERENCES 


{. BENARD, Ann. Chim. Phys. (7) 23 (1901), 62-144. 

A. PELLEW and R. V. SourHWELL, Proc. Roy. Soc. A, 176 (1940), 312-43. 
3. H. Jerrreys, Phil. Mag. (7) 2 (1926), 833-44. 

4. H. and B.S. Jerrreys, Methods of Mathematical Physics (Cambridge, 1946). 
5. K. Coanpra, Proc. Roy. Soc. A, 164 (1938), 231—42. 











som 


Joul 
It 
ang! 
the 
the 


11 
DE" 
upo 
of u 
has 
reti 
in g 
whe 
app 
grol 
doe: 
imp 
and 

A 
but 
We 
acti 
Jou 
pape 


[Quai 


50S 


6) 


ids 
the 











THE FORCES ACTING ON AN AEROFOIL OF 
APPROXIMATE JOUKOWSKI TYPE IN A 
STREAM BOUNDED BY A PLANE WALL* 
By 8. TOMOTIKA, Z. HASIMOTO, and K. URANO 

(University of Kyoto, Kyoto, Japan) 
Received 3 October 1950] 


SUMMARY 
The investigation of the manner in which the effect of the ground on the lift and 
moment of an aerofoil is modified by its thickness and camber is of great interest not 
only from the theoretical point of view but also from the practical standpoint. 


In this paper we discuss in detail the modifications introduced by the thickness 


f an aerofoil on the ground-effect upon its lift, by evaluating the forces acting on 
some aerofoils (especially approximately symmetrical aerofoils) of approximate 
Joukowski type placed in a stream bounded by an infinite plane wall. 


It is found that when an aerofoijl is placed in the vicinity of the ground at the small 


angles of attack used in practice, the lift on the aerofoil is in general increased by 
presence of the ground, but the rate of increase of the lift becomes smaller as 


the thickness of the aerofoil increases. 


1. Introduction 

DrTAILED theoretical investigations concerning the effect of the ground 
upon the lift and moment of a plane aerofoil have been carried out by one 
of us (8S. T.) and co-workers (1, 2, 3). However, since such a plane aerofoil 
has no camber or thickness, it would not be appropriate to apply those theo- 
retical results directly to the case of aerofoils used in practice, which have, 
in general, camber and thickness. Some modifications would be necessary 
when either, or both, the camber and thickness of an aerofoil becomes 
appreciable, though it may naturally be expected that the effect of the 
ground on the lift and moment of very thin aerofoils with negligible camber 
does not differ greatly from that of a plane aerofoil. Thus an interesting and 
important question arises: In what manner is the ground-effect on the lift 
and moment of an aerofoil modified by its camber and thickness ? 

A beginning has been made on this problem by Bonder (4) and Miiller (5), 
but no definite results seem to have been obtained, as far as we are aware. 
We have therefore discussed the problem in detail, by evaluating the forces 
acting on a circular are aerofoil as well as on an aerofoil of approximate 
Joukowski type with small camber, each of which is placed in a stream 


* The work described in the present paper was carried out in 1941. Publication of the 


er has been unavoidably delayed. 


[Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 3 (1951)] 
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bounded by an infinite plane rigid wall. A circular are aerofoil has been 
taken for the purpose of investigating the effect of the camber of the aero- 
foil upon the ground-effect, while an approximately symmetrical aerofoil 
of approximate Joukowski type has been considered for the purposes of 
finding what modifications, if any, are introduced in the ground-effect by 


the thickness of the aerofoil. In the present paper we describe our results 


ee Ra 





Wall 


for the latter case. Results for the former case, i.e. for the case of a circular 
are aerofoil, will be published elsewhere (6). 

Although there are several different ways of dealing with the problem, 
we have used the convenient method of images, as used in Bonder’s and 
Miiller’s papers. Thus by introducing on the lower side of the wall the 
mirror image of the aerofoil with respect to the wall (see Fig. 1) and then 
allowing the fluid to have access to the whole region, we replace the original 
problem by the two-dimensional continuous flow of an incompressible per- 
fect fluid past a system of real and image aerofoils placed in an unbounded 
stream. The required flow pattern past such a system of real and image 
aerofoils can then be deduced by suitable conformal transformations from 


the flow past a system of two equal circular cylinders. 


2. The flow past two equal circular cylinders 

Following the lines of the argument laid down by Lagally (7), we shall 
first deal with the steady two-dimensional continuous potential flow of an 
incompressible inviscid fluid past two equal circular cylinders, under the 
assumption that the line joining their centres is perpendicular to the direc: 
tion of the undisturbed flow at infinity and that there are no singularities 
in the field of flow. 
We take the plane of fluid motion as the (€,7)-plane. If Q,, Gp are 
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een the points (0,c) and (0, —c) respectively and P any point in the plane, and 
Tr if the radii from Q,, Q, to P are of lengths r,, 7, and are inclined at angles 
foil $,, 3, to the €-axis, then the complex coordinates € = €+-in of the point P 
s of are expressed in either of the following two forms: 
ul 
ra i lo . P ‘ 
and we have log log + Ut, U',). (2) 
ic / 
1 
The curves 7,/7 constant are a set of co-axal circles having Q,, Q, for 
| 2/*1 =) 1 2 


limiting points. The circles corresponding to values of r,/r, greater than 
unity lie above the €-axis and those corresponding to values less than unity 
below, while on the €-axis itself, which is the common radical axis, 73/7, = 1. 
The curves #,—#, constant are circles, or rather ares of circles passing 
through Q,, Q, and cutting the first set of circles orthogonally. On the 


right-hand side of the y-axis #,—#, is positive and on the left-hand side 





negative, while on the y-axis #,—v, 0, except on the segment Q, Qo, 


where #,—?, +a. At infinity r,/7, 1, d,—?, 0, and at Q,, Q. we 
have 7/7 © and r,/r,; = 0 respectively; also, logr,/r,, ®,—, are the 
bipolar coordinates. 
We now take a circle C, defined by r,/r, constant A (> 1) above 
the ¢-axis and another circle C, defined by r,/r, = constant = p (< 1) 
cul below, and we assume that these two circles are of the same radius a. 
Then there exists an obvious relation Au 1 between the two constants 
blen } 


and these two circles are symmetrically situated with respect to the 


ixis (Fig. 2). Further, we denote the centres of the circles by O,, O, 


I the respectively 
1 the If we denote the centre distance O, O, by 2b and if we put 
rin 
IB log A log ps e (>9@) 
le we . 
for convenience, we find 
unae 
| c bh 
image sinh o cosh oa : (3) 
a a 
s tre 
Thus when the common radius a of the two equal circles and the distance 
ipart of their centres, 2b, are given, equations (3) enable us to calculate 
the values of ¢ and a. 
1 
rT Next we conformally transform the ¢-plane on to an s-plane by the 
v ol relation e m 
ay the ° « ie 
let S hc aaa (4) 
» diret 4 uc 
larities Then the region outside the two circles C,, C, is transformed into a rectangle 
1 sides 27, 2o (Fig. 3), the circles ¢ ', Cy corresponding to the sides CD, AB 


respectively. 
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If, as shown in Fig. 2, we assume that at a great distance from the 
cylinders the fluid flows with a constant velocity U from left to right 
parallel to the €-axis, the complex velocity potential w for the flow under 


discussion takes the approximate form 





w UL, 
when ¢ is very large. 
) From the relation C c cot 3s, (5) 
| 2 
| which follows immediately from (4), it is found that € = 0 corresponds 


to s 0. Therefore, if we take account of the fact that cot 4s behaves 


as 2s when s is small, we readily find that in the neighbourhood of the 

origin (s = 0) of the s-plane the complex velocity potential w assumes the 
approximate form ; 

2cl ; 

uw a (6) 


s 





which is the complex velocity potential due to a single doublet of strength 
2cU at the origin s = 0. Since cU > 0, the axis of the doublet is in the 
positive direction of the real axis and the direction of flow in the vicinity 


of s 0 is as shown in Fig. 3. 


x : dw Pel) 
From (6) we have — ( 
ds = 


~I 
— 


which shows that the required general expression for the conjugate complex 
velocity dw/ds must have a pole of the second order at the origin s = 0. 
Moreover, it has no other poles within the rectangle, since we have assumed 
that there are no singularities in the whole field of flow in the ¢-plane. 

Our next problem is to obtain the general expression for the conjugate 
complex velocity dw/ds for the flow in the rectangle A BCD in the s-plane, 
and this can be done by the use of Weierstrass’s (-function. 

We first note that the two sides AB, CD of the rectangle A BCD, which 
correspond respectively to the circles C, C, in the ¢-plane, are both stream- 
lines of the flow: also since the homologous points on the remaining two 
sides AD, BC correspond to the same single point on the segment AD 

BC’) in the ¢-plane, the values of the conjugate complex velocity at 
the homologous points are equal to each other. Further, we recall that 
there is, as mentioned above, a single doublet of strength 2cU at s = 0 
with its axis in the positive direction of the real axis. 

These conditions can be satisfied if we take the rectangle A BCD as the 


fundamental period-parallelogram and we arrange a set of equal doublets 
| | 


I 
4 


i Strength 2cl’ with similar axes, one at the point s 0 and one at all 


the homologous points. It then follows that the conjugate complex velocity 
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for the flow due to such a system of infinitely many doublets has poles of 


the second order at s = 0 and at all the homologous points. 
Thus, introducing Weierstrass’s @-function with periods 2m, = 2z, 
2w, = 2io, we may write the conjugate complex velocity dw/ds in the form 
du ; 
2cU @(s)+k, (8) 
ds 


where k is a constant determined by the circulation round one of the circles 
(C,, say) in the ¢-plane. 

Moreover, the value of this constant can be determined if we assign the 
position of one of the stagnation points (e.g. that on the downstream side) 
on the circumference of the circle C\. 

Let the stagnation point on the downstream side on the circumference 
of the circle C, be denoted by ¢, ib+-ae~*®, Then, writing f, = €,+in, 
we have immediately 

€, = acosB, n, = b—asin§. (9) 
Further, let s, be the point in the s-plane which corresponds to the stagna- 


tion point ¢,. Then, taking (3) into account, we have from (4) 


, ib+ictae-'P ‘ ‘ sinha eos B 
8 i log - é . io—tan-! : ) (10 
ib —ic--ae-'P 1—cosh osin B 
o1 8; is +0, w.t+,, (11 
where we have put for simplicity 
sinh o cos 8 
6, tan i —al- (12 
] cosh o sin 8 
The condition that ¢ f, should be the stagnation point requires that 
- 41 1 


du: du ds ‘ 
i), G8), =* : 
dl v=t, as), 81 dl S=81 


which is also equivalent to the relation 


(| 0. (14 
\ds }. ‘“ 


since, in the present case, (ds/dZ), ,, ~ 0, as can easily be proved. 
Substituting (8) in (14), we have 
2cU ~(s,) +k (0, 
which gives k 2c —(s,), (15 
and if this value of & is substituted in (8), the conjugate complex velocit) 
for the flow in the s-plane becomes 
dw 


z 2c U{O(s)—Q(s,)}. (16 
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Further, if we combine this with (4), we can immediately obtain, by using 
the relation dw/df = (dw /ds)(ds/d@), the conjugate complex velocity dw/dl 


in the ¢-plane. 


3. Conformal transformation of two circles into two aerofoils 

Making use of an appropriate transformation formula, the region outside 
the two equal circles C,, C, in the ¢-plane can be conformally transformed 
on to the region exterior to two equal aerofoils in a z-plane, which are 
symmetrically situated with respect to the x-axis (ef. Fig. 1). A trans- 
formation formula suitable for the present problem has been suggested by 
Miiller (5), who by employing similar transformations has discussed some 
special cases of little practical interest in which aerofoils are placed at 
zero ingle O ittack. 


The suggested transformation formula is given by 
(17) 


where Cy U & 


; a m No, and & 5, no, p, y are all real constants 
which determine the shape and orientation of the aerofoils in the z-plane. 
\s will be shown later, suitable choice of these constants enables us to 
obtain a pair of equal aerofoils of approximate Joukowski type, symmetri- 
cal or unsymmetrical, which are situated symmetrically with respect to the 
axis, as shown in Fig. 1. Then if we assume that the fluid at infinity 


ows parallel to the x-axis, this axis becomes coincident with a particular 


streamline and can therefore be replaced by a plane rigid wall without 


causing any change in the field of flow. 
We first consider the transformation formula (17). Separating real and 


imaginary parts we have 


rae Ee +Fo—(n— No) tan y; 
\ So) 0)” 
' lanes sc c ti ay 
518 + Sot (9+ Ho)tan y; 
\7) No)” 
(18) 
p= cos 7 7 ) 
y se) No e-7 €,)tan y; 
\¢ Co)” 0)” 
pp- COS Y { Ps ms ) 
at) No— (+ Fo) tan ys 
( se No)” 


In order to obtain an aerofoil section with a sharp trailing edge, it is 


necessary that one of the zeros of the transformation (17) should lie on the 
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circumference of the circle C,, while the remaining zeros are enclosed within 
its circumference. 

We now assume that the aforementioned stagnation point ¢ = ¢, on the 
circumference of the circle C, is taken to be coincident with one of the zeros 
of the transformation (17) so that dz/dZ 0 at ¢ = ¢,. Then the point ¢ 
is transformed to the sharp trailing edge of the aerofoil. Since, however, 





Gl= +) 
! ib om 
! ! CEE, +eN) 
! e + & = 8+ So E 


=, 
In 
eee nner. 


J™N| 


x 





C 


Fic. 4. 
the value of the circulation round the circle C, (and consequently the circu- 
lation round the aerofoil obtained) has been so chosen that dw/d{ = 0 at 
C ¢,, the velocity at the trailing edge becomes finite and the fluid flows 


smoothly there. 


Using (17), the equation (dz/d@);_;, 0 becomes 
y2piy y2e-ty 
Lowe Do ran wae (19 
(C1 +o) (Cy Co)” 
or (C, Co)*(Cy 1 Co)” pre'M(Cy Co)” pre 'v(C, 1 i Q. (20 


It is easily found that the coefficient of ¢} on the left-hand side of (20) is 
equal to 4é,. Also, it is readily seen that the quartic equation in ¢ obtained 
by replacing f, by ¢ in (20) has two pairs of conjugate complex roots. 
Hence it will be found that the four roots may be written as ¢,, t, to Se 
and that they are situated as shown in Fig. 4. 
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hin Separating real and imaginary parts on the left-hand side of (19), we get, 


after some reduction, the following relations 


the 2 ¢ 4 2 
tan y \S1 0) (1 — No)” (€; So)" — (1+ No) | 
TOS any| 7, » i Fe 212 
1y\Sq o) (1 No) 5" \S4 €o) (1 T No) 
ver. ”¢ é ) 1 No : 1+ No | (21) 
“\31 50 (/< , > 9,2! - \o 9)91? - 
Sa So)" + (4 No)*$" iE r€&o)° + (15 No)"$" 1 
] cos — " FS “ 
( 24 | yay2\%1 £5)? —(m— no)? — 2(€, + £0)(71— n9)tan y]+ 
Pp SI 0) 1 lo) 5 
cos come " ae 
7 a 212 (€1+&o)°—(m1+ no)? + 2(61+£0m-4 No)tan y]. 
i\S] 0 Ny No 


(22) 


Thus if we assign the values of €,+& 9, 9, ,, these two equations enable 


us to determine the values of the two constants y, p, and by using these 





values of y and p we can also find the values of the remaining two roots 


In order that the transformation may be conformal everywhere in the 
region outside the aerofoil above the real axis, the circle C, must be so 
chosen that the zero-point ¢, lies within its cireumference, and this con- 
dition should be taken into account when we assign the values of the 
constants £€,;+-€, mH», 7, On the other hand, since, as we have seen, the 
other zero-point ¢, has been taken to be on the circumference of the circle 
(,, it will be'transformed into the sharp trailing edge of the aerofoil above 
the x-axis. In this case it is readily found from symmetry that the zero- 
point @, lies within the circumference of the circle C, and that the zero-point 


¢, lies on the circumference of C, and is therefore transformed to the 


arcu sharp trailing edge of the aerofoil below the a-axis. The transformation 
v is now conformal everywhere in the whole region exterior to the two 
flows aerofoils. 

We next consider the shape and orientation of the aerofoil which is ob- 
tained from the circle C, by the transformation (17). Let the chord length 
of the aerofoil be denoted by / and let the distance of the midpoint of the 
chord from the x-axis be denoted by H. Further, we denote by « the 
; geometrical angle of attack, defined as the acute angle between the chord 

” of the aerofoil and the x-axis. 
2()) 1s Now let z, Upt+iyp and z, x, +ty, be the complex coordinates of 
ained the points in the z-plane which correspond respectively to the points ¢, and 
roots cin the ¢-plane. Then z, is evidently the complex coordinate of the trailing 
Co & edge of the aerofoil. If we choose the values of the various constants so 


that the point z, becomes very nearly the leading edge of the aerofoil, the 
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distance between the two points z,; and z, may be taken approximately as If 
the chord length. Thus irrot 
] 2, — 2p V{(t,—ae)?+(y¥,—Yr)"}. (23) of n 
The distance H of the midpoint of the chord from the a-axis is given by Is - 
H = 3(¥p,+Yr), " 
, prov 
and therefore the ratio //H becomes site 
/ 24 {(@,—Xp)?+ (yz, Yr)"} 34) dens 
H YntYn % 
Lastly, the angle of attack a is given by 
y = tan (Be a (25) . 
Ln — Xz us d 
From equations (18) or from Fig. 4 it is seen that the values of /, H, and Bert 
. are completely determined if we assign the values of €,+-£ 9, 9, ,. How 
ever, it will be seen from (9) that even when the values of €,-+-&, 7, 7; are 
given, the values of the constants a and b, which determine the size and whe 
position of the circle C,, remain indeterminate. Thus it is necessary to me 
assign appropriate values to a and b. ali 
Instead of assigning appropriate values to a and b, we have adopted the 
more convenient method of assigning appropriate values to € and £ which, 
together with the prescribed values of €,-+-€, and 7,, determine easily (by 
(9)), the values of a and b. The advantage of this lies in the fact that the whe’ 
thickness of the aerofoil is intimately connected with the value of &, while this 
the value of 8 has a close relation to the camber of the aerofoil. The 
details of numerical calculations as well as the results obtained will be 
given later. D 
It may be added here that the geometrical angle of attack « is approxi- ing ¢ 
mately equal to }y, as was to be expected from the fact that the trans- we h 
formation (17) is but a slight modification of the well-known Joukowski 
transformation. 
4. The general expression for the lift acting on the aerofoil 
By combining the results obtained in the preceding sections we can 5. E 
obtain the steady uniform flow in the positive direction past two equal Fy 
aerofoils of approximate Joukowski type which are placed symmetrical to tl 
with respect to the x-axis, and since in this case the x-axis can be replaced spon 
by an infinite plane rigid wall, we can thus obtain the flow past an aerofoil s-pla 


of approximate Joukowski type placed close to a plane rigid wall. 
We now proceed to the calculation of the lift acting on the aerofoil placed 
in the vicinity of the wall. We denote the x- and y-components of the force 


acting on the aerofoil by P.. and P, respectively. 
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If, in the first place, we take into account the fact that the flow is steady, 
irrotational, and continuous, we can easily prove by the use of the principle 
of momentum that the aerofoil does not experience any resistance, i.e. P, 
is equal to zero, just as in the case of an aerofoil in an unlimited stream. 

Next, by applying the same principle of momentum, it can be readily 
proved that the lift acting on the aerofoil is equal in magnitude but oppo 
site in sign to the total fluid thrust acting on the bounding plane wall. Thus 


] 


lenoting the latter by P’,, we have 


P, yy (26) 
We shall now calculate P’. To avoid confusion with the constant p, let 


us denote the fluid pressure at any point in the fluid by p,. Then, by 
Bernoulli's theorem, we have 

Pp, = sp(U2— |v), 
where v denotes the magnitude of the fluid velocity at any point and p 
the density of the fluid concerned. Therefore, the total pressure P’, acting 


on the bounding wall is given, by 


Pp’ | pda tp | (U2—|v|?) dz, (27) 
where the integral is taken from © to © along the w-axis : and combining 
this with (26), we have 

P Lo | (U?—|v\?) da. (28) 


Denoting the lift acting on the aerofoil by L in place of P, and remember 


ing that alone the a-axis. dz lx and |v u Ch Cx Cw) ex dwidz, 
ve have, from (28 
‘ ’ w\? 
L p | ye (' ra (29) 
J | dz} J 


5. Evaluation of the lift 

For convenience in the evaluation of the above integral we transform it 
to the s-plane. Thus, remembering that the a-axis in the z-plane corre 
sponds to the real axis of the ¢-plane as well as to the real axis of the 
5 plane, we h Lveé 


] - | ly fer te) 


("-) (“*) dé. (30) 
ds dc dz} } dl r= ds}. 6 


where 6 denotes the real axis of the s plane. 
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We shall next calculate the value of the integrand. From (3), (5), and 
(16) we have 


é —asinh o cot 36, 
( 1—cos 6 
di},  asinha’ ; (31) 
lw ee 
( 4 = 2aU sinh of 9(6)—@(s,)} 
ds } 9 
Also, we have from (17) 
dz 2p" cosy 
o me f_ f(€ 1 €&,)2__ p21 Q(€4+-€,)y tany}. (32) 
ic), £ s(E T €y T ne}? : So) ” . =o) ” = 
= md dz : 
Thus, writing also : S, (33) 
dt t ¢ 
lw ds di q 
we have mt = = 2U{e(8)—(s,)}(1—cos 6)S (34) 
ds dl dz s=6 


c=é 

We now express the function (@)—(s,) in a form suitable for numerical 
calculations. Writing uw = s/2w, we have in general 

o,(s) - 


a(s) 


(o(s)—e, se 


2w, 7,(0)F,(u) 











If a complex variable v satisfies the relation 


“f)<o() <"¢h ; 


we have (8) 
l FH (0)F,(v 1 ayy) 
4n H,(v)9,(w) 


. 2n «3 » 9 9 In Dn) 
sg?" sin 2a(nv+w)—q?" sin 2nz1 
} cot zv+ } cot mw 4% 2 ( \—4 , 


ios 1 — 2q2" cos 27w-+-q*" 
n=1 
from which, by putting v = } and w = u, we obtain immediately 
lL #(0)9,(w) = gq?" sin 27u a 
71\ 5 | cot zu S (—1)" a 7. (38 
4 3,(0)7,(u) oad 1 — 2q?" cos 27ru-+-q*" 
= n=1 
In our case, w, = 7 and ws = io, and therefore, when s is equal to @, we 
have u 6/27. Thus, combining (35) with (38), we have 
: 2n oi 2 
“ q?" sin 6 , 
(o(0)—e Leot 10+2 NS ( 1)* l ; (39 
, 6 . yay | — 2q?” cos 0+-q*" 
24?" cos 0-+-q 
n=1 
The parameter q is given by q = e77'. But, since tT = wg/w, = 10/7, We 


have, using (3) b Ib 


; e | ( 
l ’ a | \ a 


} 








if w 


In 
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and Next, when s = 8; = ws+6,, we have u = 8/2w, = 0,/2w,+47. Hence, 
if we put u 6,/2w,, we have u u,+47 and 
| 1 1 1 ee 


( ‘ 
OO) 
| Vs (O)P(U, 


| F(O)F,(u, - } 


(41) 
2w, F,(0)F,(u,—4 





= 2 
51 
{)(S;) ey 
Zw, V(0)0;(Uy + $7) 


— 





In general, when a complex variable v satisfies the relation 


“f)<()<=¢} 





oz ] b (O)8. ? uw) l 
dr, (v) 4 (U 4 sin 772 
™ \ GP" sin a (2n—1)v+ 2w]+-q*" sin(2n— 1 av 
, es 2q?" leos Paw q* 2 , 
? l 
34 : , 
whence we have, by putting v 5 and replacing w by v, 
: lL #,(0)0,(v) l . q?"—1! cos 2a7v-+q4"-* 
T1Ca ] 4 % ( 1)” 1 1 1 - (42) 
tor h,(0)9.(v) } Lan 1+ 2q?"-! cos 27v-++-q*"-? 
2 3 pe | 
- By putting v = u,—} in this formula and remembering that w, = 7 and 
7 consequently wu, 6,/27, we get 
1 (OW, (u L) “ q2"-1 eos 6 an-2 
71 s\U— 2 Tw } n+ : 14 ., (43) 
. 2w, (O)Ps(U, +) 2 a l 2q"" ‘cos 0, | qi" “ 
6 “ 
ind hence. by (41). 
* qen leos A qian 2 2 
0(s,)—e, +25 (~1p* 2] 14 , (44) 
4 — 1 — 2q?"-1 cos 0, ++-q*"-? 
Subtracting this from (39) we have 
()(7) —P(s,) P+Q, (45) 
with 
= 2n 3 2 
Pp L cot 10-42 = 1)" q-" sin 6 
; “ ilo | — 2q?" cos 6+-q4" 
38 : (46) 
. ee ee 4n—2 2 
O , “ X Te q " cos 6, q we 
A. vi , — 1 — 2q?"-! cos 0,-+-q*"~ 
and this formula is conveniently used for numerical calculation of the 
" values ot OO ()(S,). 
aking (31), (33), and (45) into account, the expression (30) for the lift 
becomes. after some reduction. 
7 W 
, Ses { S " al | 2 
4( L ol’*asinho — 4(1—cos 6)(P+-Q)? —} dé. (47) 
i | cos @ S| 
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If, by using the chord length / of the aerofoil, we define the lift coefficient 


Cy, as C, = L/(}pU*l), we have 
; asinho , S <a 
( t(1— cos #)(P+Q)*—} dé. 
” l | 1—cos 0 — é) S| (ss 


Since the integrals in (47) and (48) are intractable, we have evaluated 
them numerically, and used the limiting values of the integrand when 
@ > 0, which are given by 
j_ 48 L| 


lim 4(] 0s 6)(P+Q)2 
rei cos @ \ een S| 


2 a | n qt 
l >( = cos y+ 16 S > 4@). (49) 


asinho 
4 n=1 


6. Numerical calculations 
As mentioned already, an aerofoil of approximate Joukowski type placed 


close to the ground can be obtained, in general, by transforming a circle 


in the vicinity of a plane wall by the aid of the transformation equation 
(17). But, in order to investigate the manner in which the ground-effect 
upon the lift of an aerofoil is modified by its thickness, it is convenient to 
take, if possible, a symmetrical aerofoil with no camber. Thus we have 
obtained a few approximately symmetrical aerofoils with very small cam 
ber by choosing appropriate numerical values for various constants con 
tained in the transformation equation (17). For the purpose of computing 
the effect of the camber of an aerofoil, however, we have also investigated 
some aerofoils with slightly larger camber. 

There are several ways of assigning appropriate numerical values for the 
constants in the transformation equation (17). But the most convenient 
way seems to be to assign at the outset appropriate numerical values for 
the three quantities €,+-€ , y, and y). Then the values of p and y can be 
determined by (21) and (22), and if we then obtain the second root ¢, of the 
equation dz/d¢ = 0, by making use of those values of p and y, we can deter- 
mine, by (23), (25), and (24), the values of the chord length /, the angle of 
attack «, and the ratio //H of the chord length / to the distance H from the 
midpoint of the chord to the ground respectively. 

After some trials we found that an approximately symmetrical aerofoil 
of approximate Joukowski type placed near the ground, at an angle 0! 
attack as used in practice, can be obtained by assuming &,+&) = !. 
1 1-43, and yn» = 1:57. The values of p, y, 1, « and //H in this case are 
shown in Table 1 below. 

The thickness of the aerofoil is determined mainly by the value of 
while the camber is determined by the value of 8. Since our principal 
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object lies in the discussion of the manner in which the ground-effect on 
the lift of an aerofoil is modified by its thickness, we have taken the values 
of €, and 8 as shown in the table. Aerofoils A and B are both approximately 
symmetrical aerofoils of approximate Joukowski type with very small cam- 
ber, while aerofoil C is an aerofoil of approximate Joukowski type with 
slightly larger camber. The values of €,+&, y, and 3, and consequently 


those of 7 / «, and / H as given in the table are all common to these 


three aerofoils. The values of the thickness ratio and camber for these three 


TABLE | 
/ 3°4 
5 x 5 I } 
/ B 1 ( 
1 I 
g°s I*l 
I 12° 2 
I°-11090 I*l 54 
I 1°62396 538 
137 14 
} } 17 
s saan a 
761) ( *759) ( )43) 
verofoils are shown in the table.* The values of the constants a and b as 


uculated from (9) are also given in the table. 

The aerofoils A, B, and C are shown in Fig. 5. As shown in Table 1, the 
thickness ratios of the aerofoils A and B are 0-085 and 0-137 respectively, 
while the aerofoil C has the thickness ratio of 0-144. The cambers of both 
verofoils A and B are nearly equal to 0-004, while that of the aerofoil C 
is 0-017 approximately. Since the two aerofoils A and B have cambers 
nearly equal to each other, comparison of these two aerofoils enables us 
to anticipate the manner in which the ground-effect upon the lift of an 
ierofoil is modified by its thickness. Also, comparison of two aerofoils B 
ind C with nearly equal thickness ratios enables us to anticipate the effect 
f the camber of the aerofoil upon the ground-effect of the lift. 

Evaluating the integral in the formulae (47) and (48) numerically, the 

ues of the lift L and the lift coefficient C,, of each of the three aerofoils 


lave been obt uned. To do this use has been made of the values of 6, and q 


As kness rati f aerofoil is here defined as the ratio of the maximum 
goth of the aerofoil. Also, the camber of an aerofoil is here defined 
naximu! tance f the « r line from the chord divided by the chord length. 
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which have been calculated with the aid of (12) and (40) respectively, by 


employing the values of a and 6 as shown in Table 1. 


Aerofoil A 





Wall 


Aerofoil‘B 





wall 


Aerofoil C 









wall 


Fic. 5. 
The values of the lift coefficient C,, thus calculated are given in Table 1. 
It is found that the value of C, for the aerofoil C is much greater than that 
for the aerofoil B, despite the fact that these two aerofoils have nearly equal 
thickness 1 
these two aerofoils. 








atios. This is perhaps due to difference between the cambers ol 
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te 7. The lift on the aerofoils when placed in an unlimited stream 

| In order to investigate the manner in which the lift of an aerofoil is 
changed by the presence of the ground, we have calculated the lift L, 
experienced by the preceding three aerofoils when placed in an unlimited 


stream. There are in general several methods of calculating the lift of any 





given aerofoil placed in an infinite mass of fluid, but we have employed 
the method used by Karman, Trafftz, Miller, and others.* 

The values of the lift coefficient C,, thus calculated are shown in Table 2, 
where the values of C; have been reproduced in the first row. The values 
of the ratio C;,/C,, are also given in the last row. For the sake of compari- 
son, the values of the corresponding quantities for the plane aerofoil are 


given in the last column.’ 


TABLE 2 





The results are shown in Fig. 6, where L/L, (= C,/C,,) and l/H are the 
ordinate and abscissa respectively. In this figure the solid line curve shows 
the values of L/L, for the plane aerofoil placed at an angle of attack 
.= 5 11’ 45”, which have been calculated by using the general formula 
given some years ago by one of the present writers (1). 

From Table 2 and Fig. 6 it will be seen that when an aerofoil of approxi- 
mate Joukowski type with small thickness ratio is placed in the vicinity 
of the ground at a small angle of attack as used in practice, the lift is 
increased by the presence of the ground, just as in the case of the plane 
verofoil, but the rate of increase becomes smaller as the thickness ratio 
increases. Further, by comparing the values of L/L, for the two aerofoils 
B and C’, we see that the rate of increase of the lift due to the presence of 
the ground is greatly modified by the camber of the aerofoil, the rate of 
increase being decreased as the camber increases. The latter result is in 
agreement with the result of our investigations concerning the ground- 
effect on the lift of a circular are aerofoil placed in a stream bounded by 


1 plane rigid w ill (6). 


‘able 

in that The essentia f the method are described in W. F. Durand, Aerodynamic Theory, 2 
26 SO 8 

vy equ t id of 

: ( r the plar erofoil has been calculated by the aid of the general 

i bers formula as ¢ ne of the present writers (1), assuming. as before, the value of the 


H to be 3-486 and the ang f attack a to be 5° 11’ 45”. 





X 











| 


Plane aerofoil 





e AerofoilA 
| 

x Aerofoil"B’ 
4 Aerofoil’C’ 
| 


! 

| 
| | 
| 





, 





3486 4 

















wo 


eX] 


obi 
the 
pla 
tio: 
the 
bei 


hay 








FORCES ACTING ON AN AEROFOIL 307 
Note added January 1950 


In a recent paper (10) Green has obtained the first three terms in the 


series expansion for the lift acting on the symmetrical Joukowski aerofoil 


] 

" of small thickness in the presence of a plane wall, and by comparing the 
results with those for a circular are aerofoil he has conjectured that, for 
small values of the angle of attack, the effects of camber and thickness 
work mainly in opposite directions, so that the plane aerofoil might be 
expected to give a reasonable approximation to actual cases. 

Quite recently Mr. H. Fujikawa has extended Green’s calculations and 
obtained the first five terms in the series expansion for the lift acting on 
the symmetrical Joukowski aerofoil of small thickness in the presence of a 
plane wall. Fig. 7 shows, for comparison, the results of numerical calcula- 
tions for three cases in which the thickness parameter ¢ of the aerofoil takes 
the values 0, 0-1, and 0-2 respectively, the corresponding thickness ratios 
being 0, 0-118, and 0-215 respectively. The values of the angle of attack 
have been taken to be 5° and 10°. 
6 From this figure it will readily be seen that, in accordance with the 


results obtained in the present paper, the effect on the results for a plane 
ierofoil due to thickness acts in exactly the same direction as the effect 


due to cambet 
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THE PASSAGE OF TURBULENCE THROUGH 
WIRE GAUZES 
By A. A. TOWNSEND (Emmanuel College, Cambridge) 


[Received 29 August 1950] 


SUMMARY 

Measurements have been made of the turbulent intensities existing after the pas. 
sage of isotropic turbulence through fine-mesh wire gauzes, and the results have been 
compared with the predictions of the theory of G. I. Taylor and G. K. Batchelor, who 
express the intensity reduction in terms of the pressure and refraction coefficients of 
the gauze. In agreement with the theory, the intensities of transverse components 
of the transmitted turbulence are considerably greater than that of the longitudinal 
component, but the degree of anisotropy is less than predicted. Observed reductions 
in total turbulent energy agree well with the predicted values. Upstream of the gauze, 
reductions of turbulent intensity are observed in very good agreement with theoreti- 
cal predictions. It is concluded that the Taylor—Batchelor theory represents suffi- 
ciently accurately the motion at the gauze, but that turbulent energy transfer 
occurring within the zone of influence of the gauze leads to appreciable redistribution 
of energy between the turbulent velocity components. Outside this zone, energy 
transfer is comparatively weak, and the subsequent approach to isotropy is extremely 


slow. 


Introduction 

For many years it has been usual to place woven-wire gauzes in the settling- 
lengths of wind-tunnels with the object of reducing both the spatial varia- 
tions of mean velocity due to flow around corner vanes and the intensity 
of the turbulent velocity fluctuations. The effectiveness of this procedure 
has been amply confirmed by experiment, but until recently no satisfactory 
theoretical treatment was available. L. Prandtl (1), who initiated the use 
of gauzes for the reduction of mean velocity fluctuations, gave an expression 
for the effect of gauzes on a steady stream in which the velocities along 
neighbouring stream-lines varied. In this work he took account only of the 
pressure coefficient of the gauze, but more recently A. R. Collar (2) has 
obtained a different relation. A more satisfactory theory has been de- 
veloped by G. I. Taylor and G. K. Batchelor (3), who take into account 
the change in direction of the flow when the flow incident on the gauze is 
inclined to the normal. If the angle of incidence is not too large, this effect 
can be described by a ‘refraction’ coefficient, and this coefficient and the 
pressure coefficient both appear in the final equation for the reduction 
factor. It is shown that the equations given by Prandtl and Collar are 
special cases of the general equation, representing particular hypotheses 
about the ‘refraction’ coefficient. In all theories continuity of the velocity 


[Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 3 (1951)] 
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component normal to the gauze is assumed. The theory can be extended 
to obtain expressions for the reduction of the intensity of isotropic turbu- 
lence in passing through a gauze, and Taylor and Batchelor compare their 
predictions with previously published measurements of the effect of gauzes 
on the background turbulence in wind-tunnels. These measurements were 
all made after passage of the turbulence through a large contraction of the 
tunnel, and not unexpectedly the agreement is poor. More recently G. B. 
Schubauer, W. G. Spangenberg, and P. 8. Klebanoff (4) have published the 
results of measurements made immediately behind gauzes, but these results 
are also in poor agreement with the predictions of Taylor and Batchelor. 
In particular, the result that the attenuations of the down- and cross- 
stream components of the turbulence are not significantly different is at 
variance both with the theory and with some preliminary observations 
made in this laboratory some time ago. A more comprehensive set of 
measurements is now complete, and these results will be considered in 
relation to the Taylor—Batchelor theory of gauzes and to the properties 
of isotropic turbulence. 


Notation 
Consider Cartesian coordinates with Ox in the direction of the mean 
stream and origin in the plane of the gauze. The components of the 
instantaneous turbulent velocity fluctuation are denoted by wu, v, w, with 
suffixes to indicate the values at particular parts of the field. Then, define 
U the mean velocity, 
6 the angle of incidence of the instantaneous flow, 
¢ the angle of emergence of the instantaneous flow, 
, d/6, 
p the fluid density, 
\p the pressure drop across the gauze at any point, 
K the pressure coefficient of the gauze, defined by Ap = 3pU?K. 


Mean values are taken with respect to time at points fixed in space. 


Results of the Taylor—Batchelor theory 

In their treatment Taylor and Batchelor assume that viscosity plays no 
part in the phenomena except in the immediate neighbourhood of the 
gauze, and that the action of the gauze can be described in terms of the 
two parameters K and gq, i 


. that the instantaneous pressure-drop across 
the gauze is , : ° 

Ap kK p(U+u)? 
and that the angles of incidence and emergence are related by 


d v6, 
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It is also assumed that there are no dynamical changes in the turbulence 


while it is influenced by the presence of the gauze, such as would be caused 


by the ordinary processes of turbulent transfer and decay. 


On these 


assumptions, the action of a gauze on small disturbances is linear, and the 


effect on each single Fourier component of the velocity field may be com- 


puted separately. It turns out that the reduction in a Fourier component 


of the incident turbulence 


u A, ell Ud)amy+nz 


depends only on the direction of the wave-number vector (I, 


spectral density functions Ff, G, H are defined such that 


F(l,m,n)dldmdn is the contribution to u? of components 


numbers in the range / to /+-dl,. m to m+dm. n to n+-dn. 


G(l,m,n)dldmdn is the contribution to +? of components 


numbers in the range | to /]+dl. m to m-+dm. n to n--dn. 


H(l,m,n)dldmdn is the contribution to w? of components 


numbers in the range / to 1+-dl, m to m+dm, n to n-+-dn. 


.n). If 


wave- 


wave- 


" wave 


then using the suffixes 1, 2, 3 to denote values respectively far upstream 


of the gauze, far downstream of the gauze, and immediately in front of 


the gauze, it is found that 


. _ 4 x?B2 | (]tq—ak)? 
F, = #5 or, 
45- + ( ]+ x~T K : 


2 G, H, 462+ (1lt+a+K)? 
, 1’ l \2 
and F, F, . ( X) = 
467+ (1+a+K)? 
where B? = [2(m2-+-n2)-1, 


FL (l+a—ak)?—o2(1+a4+K)? 


Knowing F,, G,. H,, it is possible to compute the turbulent intensities and 


the mean squares of the spatial derivatives, and then to find the attenua- 


tion factors of these quantities: 


9 


2 Ca yh eae 
pL U5/Uj. Vv (v3 } ws) (v7 + wy). 
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If the incident turbulence is isotropic, it follows from the condition of 


incompressibility that 
F, (m*- n*)D(1? | 9? n*), 
G, n?+-1?)D(1? + m?2-+n?), 
H, (1? +-m?)@([?+-m?-+-n?), 


and, since the attenuation factor for any particular wave-number depends 
only on the direction of the wave-number vector and not on its magnitude, 
the total attenuation factors listed above are independent of the function* ®. 
The computing of these factors is straightforward, but leads to rather 
lengthy algebraic expressions, and numerical results are given in Table 1. 
It is assumed that for A O-7, a 1-1(1+-AK)-, and that for K < 0-7, 


Y 2) K S K 


. as found by Schubauer, Spangenberg, and Klebanoff. 


TABLE |] 


‘ 

V be L bs Bs 
sc -c&&¢ o-702 
702 
5 29 "524 
$5 99 329 
SI 5f 0°229 
38 o'154 





The obvious limitation of the theory is that it neglects the dynamic 
wspects of turbulent motion and assumes that the time necessary for com- 
pletion of any of the processes of turbulent transfer and decay is very large 
compared with the time taken for a sample of turbulent fluid to pass 
through the zone of influence of the gauze. If L is the scale of the turbulent 


motion and w’ is a typical velocity, then the condition is that 


Liu D/U, 


where ) is the depth of the zone of influence. Since L and D are of the 
same order of magnitude. the condition becomes 
u'/U l. 

Ifu' and L refer to the energy-containing eddies, this condition is satisfied 
when the turbulent level is low. On the other hand, smaller eddies whose 
energy is derived from the large eddies have shorter ‘reaction times’, and 
may adjust themselves to the changed conformation of the turbulence suffi- 
ciently rapidly to make the overall effect of the gauze on them appreciably 


The paper of Taylor and Batchelor contains a contrary statement, but after publication 
the error wa ticed by the autho Che factors computed are therefore of general validity 
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different from the theoretical prediction. These simple considerations 
lead to two necessary conditions for the validity of the theory, i.e. 

(a) w'/U <1, 

(b) the eddies concerned shall not be too small compared with the 


energy-containing eddies. 


Experimental arrangements 

The measurements were made in the small wind-tunnel in the Cavendish 
Laboratory, which has a sufficiently low level of free-stream turbulence for 
this type of work. The greater part of the measurements were of the 
passage of approximately isotropic turbulence through gauzes of various 
solidities placed 79 cm. down-stream from the turbulence-producing grid. 
Apart from a small number of special measurements, the air-speeds were 
480 cm. sec.~!, 620 cm. sec.-!, and 905 em. sec.-!, and the approximately 
isotropic turbulence was produced by passage of the air-stream through 
bi-plane grids of periodic spacing 2-54 em. and 5-08 em. and of spacing 
diameter ratio 5-33. 

The gauzes are of woven wire, fixed in light-alloy frames under sufficient 
tension to keep them reasonably flat. In the absence of incident turbu- 
lence, the flow down-stream of the gauzes changed from steady laminar 
flow at low wind-speeds to a three-dimensional periodic flow at inter- 
mediate speeds and finally to fully turbulent flow at higher speeds. Behind 
the higher resistance gauzes some spatial variation of mean velocity could 
be detected, approximately +3 per cent. from the mean for the highest 
resistance gauze. The effect of these variations on the transmitted turbu- 
lence is discussed below. The properties of the gauzes are given in Table IL. 

TABLE I] 


Wire eer he 
Gauze Mesh liameter | Solidity v l 480 | U = 620| l gos 
mm. mm. 
A 0°635 234 0°60 26 3°4 3°2 31 
B 0°508 O°152 o's 48 2°31 2°08 1°87 
C 0°705 O°173 0°43 50 1°44 1°32 1°21 
D 1°59 0254 0°204 53 0582 0°535 506 


Measurements were made of the mean squares of the down-stream 
velocity component, u®, and its derivative (éu/éx)? in front of the gauze, 
and behind the gauze the quantities w?, (@w/éa)?, v2, (@v/éa)? were measured. 
A few determinations of the spectrum functions of u and v were also made. 


* The critical Reynolds number is based on the wire diameter d and the critical ail 


velocity, U., at which the flow behind the gauze first becomes periodic. For an infinite 


cylinder, U, d/v 40 (Kovasznay 5). 
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itions Standard techniques of hot-wire anemometry were used for all these 
measurements. 
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iedaie Experimental results 
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ai turbulence changes from smooth steady flow to a regular periodic wake 
llat . 2 


motion and then to fully turbulent flow as the wind-speed is increased. In 
ll arrangements the scale of the wake motion behind the gauze was much 


smaller than the scale of the incident turbulence, and it seems likelv that 
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the two motions are nearly independent, i.e. the intensity due to the 


transmitted turbulence can be obtained by subtracting the turbulent 


intensity due to the gauze alone from the total observed intensity. In this 
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way curves of the type shown in Figs. 1-4 were obtained, showing the 


variations of turbulent intensity on both sides of the gauze. It is noticeable 


that the two gauzes of high resistance produce far down-stream a lar 


cP 


increase in the intensity of the down-stream velocity fluctuation, but only 


when an appreciable level of turbulence is incident on the gauze. This 
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phenomenon was observed by Schubauer, Spangenberg, and Klebanoff, 
and attributed to convection by the turbulent motion of the mean velocity 
pattern produced by local irregularities in the gauze. This explanation is 
consistent with the relative absence of any similar effect on the cross- 
stream component, but there is some evidence that the effect is an initial 
stage in the laminar instability of the mean velocity pattern, and that the 
turbulence provides an initial perturbation that is amplified. The two 


lescriptions are identical in the initial stages of development, but the 
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Fic. 5. Attenuation factors for the longitudinal and transverse 
velocity components. 


ecasional appearance of a small increase in the turbulent dissipation sug- 
gests that the instability hypothesis is more appropriate. For the practical 
problem of determining the reduction in intensity, the phenomenon can 
e disregarded if measurements are taken within 25 em. of the gauze. 

Fj 


ongitudinal intensity, and v, the reduction factor for the transverse inten- 


g. 5 shows the collected results for up, the reduction factor for the 


sity. The actual points represent the ratio of the intensity at a point 
behmd the gauze to the intensity at the same point with the gauze 
removed, plotted against A’, the differential pressure coefficient of the 
gauze. K’ is defined by 


( 


(Ap UK’ 
7 Ap) pUK 


ind is related to K by 


A K{14 | o(log a 
2 d(log U) 
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This differential pressure coefficient is the one considered in the theory, 
Generally the values of v are less than those predicted by the theory and 
the values of » are greater. If values of 4(u+2v), the total attenuation 
of turbulent intensity, are compared with theoretical predictions the agree- 
ment is good for large values of K’ but poor at small values (Fig. 6), where 
the observed factors are significantly less than those predicted. 

The turbulence behind the gauze is distinctly anisotropic, but the re- 
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Fic. 6. Attenuation factors for the total energy. 


covery of isotropy is very slow in the region accessible to observation (see 
Figs. 3-4). On the other hand, there is a rapid approach to local isotropy 
in the sense that, closely behind the gauze, 


ou\? 1 {dv\? 

(-*) : sla) ; 
within experimental error. This became evident in attempts to measure 
the attenuation factors for these quantities, ie. uw’ and v’. The actual 
values of these attenuation factors change rapidly with distance of the 
point of observation from the gauze, but at any particular point the values 
of »’ and v’ were always equal within experimental error. Since in isotropic 


ou\? 1 (ev\2 
ox 2\ox} ’ 


this must equally be true in the anisotropic turbulence behind the gauze. 


turbulence 


Some trouble was taken to obtain measurements of turbulent intensities 
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on the up-stream side of the gauze, for it seems likely that the idealized 
theory should correspond most closely with practice in a region where the 
‘xact nature of the wake flow is unimportant. A hot-wire was mounted 
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on a holder that projected from the side of the wind-tunnel and could be 
rotated about an axis in the tunnel wall. By careful adjustment it was 


possible te 


measure u° and (Cu cx)* at points within 8 em. of the gauze. 
Both these quantities vary very rapidly close to the gauze, and much 
depends on the nerves of the experimenter who adjusts a delicate hot-wire 
to come within 1 mm. of the gauze with the knowledge that contact will 


break or damage the wire and so invalidate previous measurements. For 
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this reason measurements of 43 can be no more than upper limits to the T 
true value, and only measurements of yg are given in Fig. 8. These values pred 
are very close to the theoretical values, and while the results for p3 are prisi 
very scattered their range is such that it is considered possible that experi- ner 
ments on a larger scale might verify the theory quite closely even in this om 


respect. These results may be summarized. 

(a) The reduction of turbulent intensity on the up-stream side of the gauze 
is described by the Taylor—Batchelor theory with adequate accuracy. 

(b) The reduction of total turbulent intensity in passing through a gauze 
is consistent with the theory except for small values of A’, when the 
observed reduction is greater than is predicted. 

(c) As predicted by the theory, the turbulence emerging behind the 
gauze is anisotropic, with the intensities of the lateral components 
exceeding that of the longitudinal component. The degree of aniso 
tropy is less than the theory predicts. 


(d) The recovery of isotropy down-stream of the gauze is very slow, but 





local isotropy is valid within the range of observation in the sense 
that 


Discussion 

The experimental results outlined above bear on two distinct problems. | 
The first and more practical problem is that of describing the passage of 
turbulence through a fine-mesh gauze, and it appears that the Taylor 
Batchelor theory gives an accurate description of the behaviour up-stream 
of the gauze. Down-stream of the gauze the correspondence of theory and _ beco 


observation is not so close, and two sources of divergence may be dis- — be s¢ 
tinguished. There is the excessive reduction of total turbulent intensity — flow 
behind the gauzes of low resistance, and the tendency for the turbulence — 7? 
to be less anisotropic than the theory demands. of is 


The excessive reduction of total turbulent intensity is unexpected. but _ by tl 
it is possible that the results for the gauze of lowest resistance are in- 
fluenced by the comparatively large mesh of this gauze. When flow behind 
the gauze is turbulent (as was true for gauze D at all wind-speeds), the 
possibility exists that gauze turbulence may act as an eddy viscosity to | but t 
increase temporarily the energy dissipation of the transmitted turbulence. | to re 
and so to increase the apparent reduction factor. Alternatively the passage | cones 
through the gauze may distort the structure of the turbulence in such @ | spect 
way that the rate of energy dissipation becomes abnormally high for a | at ze 
short time. The first explanation is probably correct, at least for the gauze | ratio 


of lowest resistance. initia 























THE PASSAGE OF TURBULENCE THROUGH WIRE GAUZES = 319 
the The tendency for the transmitted turbulence to be less anisotropic than 
‘Ines redicted is very natural, and, qualitatively, not surprising. What is sur- 

an rising is that, after a short period of comparatively rapid interchange of 
er] nergv between the velocity components. the ratio of the intensities be 
p l : . 
this omes very nearly constant, and the once powerful tendency to isotropy 
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re 
v an ecomes a weak drift barely detectable (Figs. 3-4). Some explanation may 
e dis be sought in the comparatively low Reynolds number of the turbulent 
nsit flow down-stream of the gauze, and the initial rapid decrease of the ratio 
ilenct ~u® could be caused by adjustments of the smaller eddies in the direction 
f isotropy. That these adjustments occur, and occur rapidly, is shown 
1. but vy the validity of the local isotropy condition 
re 
; i" 1 (cev\* 
yeh > 
’ 2\oxr] 
sty but the measurements behind the gauzes of highest resistance are difficult 
ilence to reconcile with the notion that the initial approach to isotropy is purely 
assagt neerned with the smaller eddies. For example, measurements of the 
sucl spectrum of w in Fig. 9 show that the magnitude of the spectrum function 
. for t zero wave-number is reduced in the ratio 0-090. The theory predicts a 
gauze itio of 0-022. so even the large eddies determining this quantity have an 


) 


itial rapid approach to isotropy. It is possible that flow behind the gauze 
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is not described by the theory at all, but this is improbable in view of the 
accurate confirmation of the predictions of total energy reduction. A more 
likely hypothesis is that, while the turbulence is in the zone of influence 
of the gauze, its configuration is such that interchange of energy between 
velocity components is very rapid. On this view, the initial rapid readjust- 
ment occurs only in the zone of influence of the gauze, and the subsequent 
slow approach to isotropy of the energy-containing eddies is typical of free 
turbulence. 

This leads to the second problem, the nature of the processes which tend 
to establish isotropy in a field of homogeneous but anisotropic turbulence. 
It may be assumed that the turbulence behind a gauze is axisymmetric, 


which is the next simplest form of turbulence to isotropic, and details of 


the inertial transfer of energy between the velocity components might be 
obtained by measuring the approach to isotropy. The slowness of this 
approach makes it impossible to do anything quantitative by this method, 
but this slowness of recovery of isotropy is in general agreement with results 
of recent detailed studies of the rate of inertial redistribution of energy in 
the spectrum of isotropic turbulence (Stewart 6, Stewart and Townsend 
7), and of the energy balance in turbulent wake flow (Townsend 8). All 
this work tends to show that in free turbulent flows, i.e. flows not in any 
way confined by rigid boundaries, the customary image of feverish activity 
and instability of the energy-containing eddies is very far from the truth. 
and that these eddies are in fact highly permanent structures communi- 
cating energy preferentially to smaller eddies which contribute little to the 
total energy. This is clearly consistent with the present observations which 
show that only the smaller eddies have a reasonably strong tendency to 
isotropy. It would be very interesting to repeat and extend the work of 
MacPhail (9), who distorted a field of isotropic turbulence and observed 
the recovery of isotropy. Such an experiment avoids the difficulty of the 
present experiments, that an appreciable degree of anisotropy involves 3 
large reduction in total turbulent energy. 
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THE ROUND LAMINAR JET 


By H. B. SQUIRE 
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SUMMARY 
An exact solution of the equations of viscous fluid flow for axially-symmetric 
motion is derived. It is shown that this corresponds to the round laminar jet, or, 
alternatively, to the flow produced by the application of a force at a point in a 
viscous fluid. Some examples of the calculated streamlines are given. The effect 


f a source of heat at the point of application of the force is also considered. 


1. Introduction 

THE present paper discusses an exact solution of the Navier-Stokes equa- 
tions for the flow of a viscous incompressible fluid for the case of axial 
symmetry. A special form of the stream function is assumed which corre- 
sponds to the component velocities being inversely proportional to the 
distance from the origin. The calculated flow turns out to be the exact 
solution for the flow in a round laminar jet emerging from an orifice at or 
near the origin; this flow has already been studied by Schlichting (ref. 1, 
section 57) using the approximations of boundary layer theory. Itis further 
shown that the flow may be regarded as that resulting from the application 
of a force applied at a point in a viscous fluid which is at rest at infinity. 
If, in addition to the force, a source of heat is introduced at the same point 
we obtain an idealized representation of a heated laminar jet. 

From a mathematical viewpoint, the flow discussed here is the axially- 
symmetric analogue of the two-dimensional flow between non-parallel 
plane walls (ref. 1, section 42, or ref. 2). In both cases the velocities are 
inversely proportional to the distance from the origin and this condition 
could be regarded as the starting point of the analysis for each case. 

The relation between laminar and turbulent jets is briefly considered 


in the appendix. 


2. The equations of motion 
We take spherical polar coordinates (r, 6,4) with @ measured from the 
axis of the jet (Fig. 1). The velocity components are (u,v, 0) measured in 
the directions of (r,@,¢4) respectively, and all quantities are independent 
of d. The equation of continuity is (1, section 41) 
S Pia ee a 
r= or rsin 6 c@ 
(Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 3 (1951)] 
5092.15 Y 
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and the equations of motion are 


cu. vou wv 
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where p is the pressure, 


V2 


ine of Round Laminar Jet 
F/pv? = 34 76 
Fic. 1. 


p is the density, v is the kinematic viscosity, and 


Og. « l eae c 
= — $90. fa}. —_ sin @—}, 
ror\ oO r?sin 6 c@ o6 


We shall assume that the stream function % has the form 


so that 


us vr f (8) 


I Cus v - 
=: — f'(9) 
r?sin@c@~ rsin@ t 
l Cus V | : 
rsin @ cr rsin @° 


Thus the velocity components are inversely proportional to r and hence 
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r*sin @ c6@ cb r*\o0" — cé 
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and, from the equation ol continuity, 


Cv 
7 !weoté 0. 
oo 


‘J 
cot @— v cosec?4 0. 
ob c= ov 








- Using these relations the equations of motion become 
u-ty- vou l cp Vv C , CUu\ , 
4 Po — sin 6 . (3) 
r r 6 por r*siné cé c6 
al v Ov lep,v cu 
‘ Po (4) 
ro6 pr cb6 hr" 06 
From (4) we obtain the pressure by integration as 
- p Po yp Vu Cy ” 
a oe eee (5) 
p 2 r+ 
| where ) is the pressure at infinity, at which wu and v vanish, c, is a constant 
of integration, and we exclude any other terms which depend on r only by 
consideration of (3). Hence, from (5) and (2), 
l Op = 2vu 2c, 
A Or > 2 73 
and (3) becomes 
u> veu 1 8 cu 2c, ; 
2u +-— sin @ +—. (6) 
} r oO ia sin 6 06 fal] © 
. 
Before substituting into (6) we change the variable to » = cos@, so that 
sin@ = ,/(1—p?) and 
1 d d 
sin 6 dé du’ 
when equations (2) become 
l Vv ; AL 
u f (we) v aa tat ) —. (7) 
2 7 r vy (I pL") 
Substituting in (6) we obtain 
og a aia d i ‘ 
hent P(e) P+ POSH) = 2f'(H) + 5-[— 2?) f"(n)]— 2. (8) 
du 
Integration gives 
ff’ 2f+-(1—p?)f" —2c, p—Ces, (9) 


where c, is a second constant of integration. Integrating again we get 
f tuf+2(1—p?)f’—2(c, p74 Cop+Cs), (10) 


where ¢c, is a third constant of integration. 
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We now suppose that (10) has a solution of the form 
ti x(1+-p)+B(1—p), 
where « and f are arbitrary constants which depend on ¢,, ¢,, and cy. (The 
existence of a solution of this type requires that one linear relation between 
€,, C9, and cy shali be satisfied.) Then a more general solution of (10) can 


be derived by standard methods to be 


pL 
a fis p)e “(1 pe 
1+ p)+A(1—p)-+2(1—p2 = d 
f = a(1+p)+B0.—p)+2( welt 1 / a a wu)? a. 
1 


where a is a further constant of integration. For the remainder of the 
present paper we shall restrict consideration to the case « — 8 = 0 for 
which this solution becomest 

f 2(1—p?*) 2 sin?0 


as 11) 
a+-l—p a+-l1—eos@ \ 


3. Description of the flow 
The streamlines are given by 
ub vr f (w) constant, 


where f(y) is given by (11). They all have the same shape for a given 
value of a, and are shown for a 1, 10-1, and 10-7 in Figs. 1-3. 

The stream tubes have a minimum area or throat at the station which 
is given by p = cos@ = (1+-a)-!, and this direction is indicated in Figs. 1-3. 

We can interpret the solution as a jet issuing from a nozzle which 
coincides with one of the streamlines as far as its throat; it must be 
assumed, however, that a special frictional boundary condition is satisfied 
on the walls of this nozzle. In the central part of the jet the flow is very 
similar to Schlichting’s approximate solution which can be obtained by 
assuming the pressure to be constant, by neglecting the first term on the 
right-hand side of equations (6), (8), (9), and by supposing that (l1—p) <1. 
In place of (11) we should obtain 


f(y) = — —, (11a) 


valid for (l—p) < l,a l. 

There is, however, a more fundamental interpretation of the solution. 
This is obtained by supposing that the flow is generated by the continuous 
application of a force F at the origin acting in the direction of the resulting 
jet. The magnitude of this force can be evaluated by integrating the rate 


t This also corresponds to c, = cy = ¢ 0. 
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of transfer of momentum across a sphere 7 = const. due to the pressure 
and kinetic momenta and to the viscous stress. We obtain the result 


F 32(a-+1) 


1 @Q 2 ¢ ‘ 
——,=5 — +. 8(1-+a)—4(1-+-a)* log| (a--2)/a]. (12) 
277 pv 3a(at-2) 
The most important contribution for a | comes from the jet momentum. 


The following numerical values were obtained: 


F'/ pv? 3°476X I0 314 107 3°282 x 10° 
[t is interesting to note the importance of the non-dimensional quantity 
| force/pv?| in the present problem. 


4. Temperature field of a heat source at the origin 


Suppose that, in addition to the force F, we have a heat source of 


strength Q at the origin, such that its effect on the velocity and density 
distributions may be neglected except very near the origin. This system 
can be realized experimentally by a hot jet with a small orifice at or nea 
the origin. Then the temperature at any point in the field satisfies the 
equation (1, chap. 14) 

oT vot se 

— -| KV2T, (13 

or r oo 
where 7’ represents the excess temperature above the temperature of the 
fluid at a large distance and « is the thermometric conductivity. We put 


7 l 
f | q(@), (14 
: 


so that, with V2 as in section 2, 


ome l of. oT 
V27 — g (sin ' 
r=sin 6 06 fale 


Changing the variable to » = cos@ and substituting for the velocity com 
ponents from (7), equation (13) becomes 
I 1 


’ , ld i 
P(e )g(e) HF (Hg (He) [(1—p2)q’(p)], 
oO du 
where o = v/x is the Prandtl number of the fluid. On integration we obtain 


ofg = (l—p?*)q’, 
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sure the constant of integration vanishing since f = 0 on the axis of the jet 
| 1. Substituting for f from (11) we get 
] 20 
) > 
2 q a+1l—p 
tun ar 5 i — - T ae comeneieen 
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Fic. 4 
the solution of which is 

+hg a 20 - 
of tl “ Al ; (15) 
e put a+l1—yp 


where A is a constant which depends on the source strength. The source 
strength Q is determined from the heat flux across a sphere of radius r, the 
quantity of he: 


t crossing unit area in unit time being pc, v, 7T', where v,, is 


n 


the velocity normal to the area and ¢, is the specific heat of the fluid at 
p 


constant pressure. In this way we obtain the formula 


1 
0) — 
= A | f'(u)g(p) dp, 
2a pC, 1 J 
ind, with f und g(x) given by (11) and (15), 


@— = al(St*\a_-(2y"\-(a\p-(a) 4]: 




















328 H. B. SQUIRE 
approximately, and then, from (14) and (15), Th 
the 
sida 2o0+-1)Q a 20 r 
| i ; : ° (16) * 
Spc, vr \a- l—p IS I 
lan 


Fig. 4 gives the temperature contours for a 0-01, corresponding to the 





velocity contours of the jet shown in Fig. 3, taking o = 0-72 (the value | 


for air), and ra. 
2o0+1)0 2. 

( “ M4 A = 10. 3. 
S77 pe pv ] ; 


We may note that for a 
given by (7) and (8), is 


— 
~ 
= 


1, the radial velocity distribution, 4. 


Comparing this with (16) we see that the temperature and radial velocity 





distributions are similar within the jet for o = 1. : 
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APPENDIX 
The Relation between Laminar and Turbulent Jets 


Though very different in detail, the flow in a turbulent jet bears a general resem- 
blance to the flow in a laminar jet. As will be seen from the analysis below, th 
turbulent jet is an example of turbulent flow for which the eddy viscosity may be 
treated as a constant throughout the jet, at least as far as the main features ar 
concerned. 

We first determine the value of a which corresponds to the turbulent jet, from th 
observed rate of spread of the jet. Experiments (3, 4) show that the stations at 
which the velocity is one-half the velocity on the axis at the same cross-section lit 
on a cone of semi-angle 5°. From the theory of the laminar jet given above we obtain 


this observed rate of spread if we take a 0-0092. Hence the flow in the jet, and 
more particularly the flow outside the jet, will resemble that shown in Fig. 3 for 
which a 0-01. If € is the eddy viscosity of the turbulent jet then, from (12) witl 
a 0-0092, 
F ee, 
= 476. 
27rpe* 


Also if the jet is formed by a uniform stream of velocity U, coming out of an orifice 
of diameter d the applied force F is equal to the jet thrust so that 


F ] mpl j2 2, 


and hence € 0-0162U,d. 











H 


weit 


iV1sl0 
nme! 


of the 





THE ROUND LAMINAR JET 329 








This equation gives the value of the eddy viscosity of the turbulent jet in terms of 
the exit conditions 

The mean velocities for the turbulent jet can be calculated from equation (7) if v 
is replaced by e. It is shown in (3) that the velocity distribution calculated for a 


laminar jet n good agreement with measurements made with a turbulent jet. 
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THE EFFECT OF FRICTION ON SURGES 
IN LONG PIPE-LINES 


By A. M. BINNIE (Trinity College, Cambridge) 
[Received 16 January 1951] 


SUMMARY 

A pressure fluctuation occurs in a pipe-line when the velocity of the liquid in it is 
altered by means of a valve. A method of calculating these surges in a gravity pipe- 
line is described for the two cases when the exit velocity is reduced to zero (i) instan- 
taneously and (ii) at a uniform rate. Allowance is made for the effect of friction, 
which is assumed to be proportional to the velocity, but the method may be used 
to provide an approximate result when the friction loss follows some other relation, 
Dimensionless quantities are employed, and the working is expressed in terms of a 
single characteristic number k, which depends upon the properties of the pipe-line 
and of the liquid in it. The method holds good unless k is too large. For the two 
cases, numerical results are given for a small and a large value of k, and they ar 
compared with those derived by the graphical method of solution. 


1. Introduction 

WHEN a valve controlling the flow of liquid in a pipe-line is moved, a 
pressure surge occurs which may have serious consequences. Hitherto the 
prediction of these surges has commonly proceeded on the assumption that 
friction effects are negligible. If the pipe-line is long, this supposition is 
doubtful, but the matter has received no great attention, and few reliable 
experimental results have been published. However, a report issued by 
the Joint Surge Conference (1) now makes available a full account of the 
pressure fluctuations observed in very long oil pipe-lines by means of 
modern electrical gauges linked together by wireless. Of particular interest 
are the results obtained with gravity pipe-lines when the control valve at 
exit was shut either suddenly or in such a way that the velocity near the 
valve was reduced at a uniform rate to zero in a known time. A graphical 
analysis is given, which (as explained in reference 2) is based on the work of 
Allievi, Bergeron, and Schnyder: this allows for friction, but supposes that 
it is concentrated at a finite number of points. In the present paper this 
limitation is removed, and a way of calculating the surges is described, 
the friction being distributed uniformly along the pipe-line. The labou 
involved is trivial compared with that of carrying out the experiments. 
In general, the method will not hold good until the motion ceases, but it 
will usually be applicable at least to the first surge in which, of course, the 
greatest pressure is developed. 


(Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 3 (1951)] 
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The same idealizing assumptions are made as in the graphical method. 
The liquid, the pipe, and its surroundings are taken to possess perfect and 
linear elasticity so that the velocity of a pressure wave is constant. The 
reflections at the open and the closed ends of the pipe-line are also perfect, 
and the motion is regarded as one-dimensional, the velocity over a cross- 
section being uniform. The only loss of energy which is taken into account 
is the same as that shown up in steady flow by a fall in pressure along the 
pipe-line. The particle velocity is supposed small; therefore the product 
terms in the equation of motion may be neglected. The surges under con- 
sideration are analytically similar to those in an electric cable, and the two 
problems mentioned above are solved by operational methods developed 


from Heaviside’s classical work. 
2. General theory of the waves 


To avoid complications arising from the varying gradients of actual pipe- 


lines, we consider the arrangement depicted in Fig. 1. The uniform pipe- 





I 











a Z a 
re a “7 u 7 ea 
a = a 





Fic. 1. Arrangement of pipe-line. 


line, of length / and internal radius a, is assumed to be horizontal; it is 
supplied from a large reservoir, and the flow is controlled by means of a 
valve at outlet. The displacement and mean velocity at a cross-section 
distant from the valve are denoted by € and u, which are connected with 


time ¢ by the relation 


0& 

u mi 2.4 

ot \ 

We also have that S of (2.2) 
Cx 

p = po(1+s), (2.3) 

P C*p9 8, (2.4) 


where s is the condensation, p is the density, c is the velocity of sound, 
ind P is the pressure measured from a zero where the density is po. 
lhe dynamical equation for an element of unit length (Fig. 1) is 
ou : 1oP 7 c? ¢ 
x Rane 


+] 
ct p Cx Po Cx 
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where X is the extraneous force per unit of mass. No information is 
available by means of which X can be determined when w is varying, and 
we must be content to employ the value which obtains under steady con- 
ditions. In steady laminar flow the tangential stress at the wall is 4yu/a, 
where p is the viscosity, and it acts over an area 27a. The mass of the 
element is pza*, hence 
Ku, say, (2.6) 
where v is the kinematic viscosity. When the initial motion is turbulent, 
the friction loss varies approximately as u?, and even with the simple index 
2 the ensuing equations become too complicated to be of service. To avoid 
this difficulty, the factor in (2.6), by which wu is multiplied, will be adjusted 
to give the correct loss of head at the pipe exit under the steady conditions 
which prevail before the valve is moved. This procedure seems the best 
that can be adopted, and it is reasonable in view of the uncertainty men- 
tioned above. We then find with the aid of (2.3) that (2.5) becomes 

a Ku—2 & 

ot Ox 


(2.7) 


Since numerical work will be required later, it is desirable to remove 
dimensional factors from the equations. This can be effected by writing 


, , , l , l r 9 ' 
2 = fe’, e = &, “= cH’, t t', k=-K, P= Cet’, 
Cc Cc 

2.8) 

in (2.1), (2.2), (2.4), and (2.7), which then become 
ui = =, (2.9) 

Ot 
ee ae (2.10) 
Cae’ 
r S (2.11) 
ead — —ky’ = (2.12) 
ot Cx 


Only one number k remains, and this may be termed the characteristic 
number of the pipe-line and the liquid in it. Hereafter the dashes will 
be suppressed, and all quantities, including the boundary conditions which 
will be introduced, will carry non-dimensional significance. 

To determine the consequence of the closure of the valve, it is necessary 
to form the differential equation for €. This is done by substituting (2.9) 
and (2.10) in (2.12) with the result 
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3. Instantaneous closure of the valve 


ind We suppose that before the valve is shut the steady velocity throughout 
on the pipe-line is w U’, the variations in velocity due to the changing 
tla. condensation along the pipe-line being negligible. The condensation at the 
the inlet to the pipe-line is denoted by S, which remains constant throughout 


all the ensuing changes in wu. Then from (2.10) and (2.12) 


2.6) ‘ Sa+kUa—tkU2?, (3.1) 
ent, an expression which is arranged to give s = S at x l and é= 0 at 
dex z= 0. This equation and 
void = U (3.2) 
sted ot 
lons are the conditions to be satisfied at f 0. 
best If the Laplace transformation be defined by 
nen 

ij( p | ely (t) dt, (3.3) 

‘ 0 


the subsidiary equation corresponding to (2.13) and subject to (3.1) and 


nove (29 


(9.2) 1S 
os ds , ‘ , 
iting .— gf U—(p+k)(—Sa+kUx—tkvU2?), (3.4) 
ax" - 
Ws | ‘ 2 k D 
where q p(p+k). (3.5) 
(2.8 The solution of (3.4) is 
U — Sa+kUx—ikU 2? 
é A sinh ga+ B cosh qa ——, 4 = : (3.6) 
2.9 p- Pp 
A and B being constants to be determined from the boundary conditions. 
2.10 \fter the valve is shut, the displacement there remains zero and the con- 
| 
2.11 densation at inlet retains its value S. Hence for ¢ > 0, 
QO atz (), (3.7) 
Ae 
. dé S 
= ind, from (2.10) at x ‘. (3.8) 
-TISTL da Pp 
S will » wate ° . . 
hict On inserting these conditions into (3.6) we find that 
whic 
dé Uqsinhg(ii—ax) , S—kU+kU2 (3.9) 
gsa\»y da p~ cosh q Pp ; a 
y (2.9 Now nl 
in a vr) , 1x(] e-2a(1 2))(] 1_¢ 2a) 1 
cosh q 
| ) 13 LZ 


ed e-Ke-xr y Q2+.2r)_| ¢ q4—zx)_} — (3.10) 
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and from Carslaw and Jaeger’s analysis (3) of electric transmission lines 
it can be deduced that the inverse transform of 


© gs 
p* 
: 

is [et Ak(P2—a2) +k | eM Ak(72—22))} dr]H(t—a), (3.11) 


in which Heaviside’s unit function 
H(t—T)=0 (<T))\ 
o (3.12) 
1 («>T)) 


and J, denotes the Bessel function of the first kind. Thus we find from 
(3.9) and (3.10) that 


i s= u([: Akt T $3 k(t? x?) Lk | ¢ WrIy(Ah(o2— 22) dr | HU x) 
JeMUy(Ne(®@ 2a) RK [ eH (Ahr? (2—a))9 dr |» 


H(t—(2 2) eb) +(SHRU FEU), (3.13 


It will be noticed that the final term represents the pressure which prevails 
before the valve is shut. 

A check can be applied to (3.13) by ascertaining the value of P when 
to. Carslaw and Jaeger (3) proved that, when ¢ — «, (3.11) tends to 


»\1 1 
kV) arte, (3.14 
R 


where R is a constant which in the present problem tends to zero. Writing 
for brevity « = (kR)!, we see that the Bessel terms in (3.13) reduce to 


(e-% e—X(2—2x) e—K24+2)_1 p—a4—z)_| eh (3.15 


The alternate terms form two geometric series; and when these are summed 


(3.15) becomes 


y I: e— we e—(2—2) | : ear — e—X(2—22) ; 
Uh poi Uk LU ; . (3.16 
i ‘ Le 


When a —- 0, (3.16) tends to 
Uk(1—2), (3.17 


therefore, from (3.13), P — S as it should. 
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To avoid the troublesome integrals in (3.13) we can expand the Bessel 


functions and take » ot 


[,(z) eat se (3.18) 
4 64 

using only the first three terms. The numerical results set out in section 5 
make it clear that 3k may be expected to be substantially less than unity, 
and normally we shall not need to consider values of ¢ exceeding 3 or 4, 
thus the argument of the Bessel function is unlikely to be much in excess 


of unity. Correct to four fig 


ures, the value of J,(1) obtained from (3.18) is 
1-266, which is the same as that given in tables, hence the accuracy of the 
approximation (3.18) should be ample for practical purposes. It is simpler 


to insert (3.18) into both parts of (3.11), rather than in the integral alone; 


and when this is done, the result is 


ly kar k ) 1 J29-21 p t 3 h:2y2 ! 1 keAg4 | 


t 8 ( J 4 4 : 1024 


It ( -h:2a?) j-242( gisk*x) | ake | oA H(t x). (3.19) 


The exponential terms can, if desired, be expanded to yield a formula 
expressed solely in powers of kx and kt, but for numerical work it is easier 
to employ (3.19) as it stands, in the later terms x being replaced by 2—z, 
2+2, 4—2z, etc. 

A crude approximation, causing an immediate error of rather more than 
| per cent. when z is unity, is to employ merely the first two terms in (3.18). 


The expression corresponding to (3.19) is then 


(3+ dhx}—e—- {2 — 2 ka? + bht+4 0} A (t—z2). (3.20) 


16 
The maximum pressure developed in the pipe-line may readily be ob- 
tained. The most heavily loaded part is at the valve, and we see from 
3.13) and (3.19) that for x = 0,0 <t < 2, 


P Uis—, kt(11 If kt? L8 | vibak't4)} L(S—kU). (3.21) 
When the exponential term is expanded, (3.21) becomes 


P = U{1+}ht—4,? + YP —aigk*t4}+ (S—kU). (3.22) 


U6 
Thus, as we expect, the initial pressure-rise above the working pressure 
is U. The pressure continues to increase until ¢ 2 when 


P= UN1—1}2438— 345, (3.23) 


1 


therefore the maximum pressure-rise above the static pressure is diminished 
by the fi iction ee , P ‘ 
) k2— LIS Bok, (3.24) 
compared with the rise calculated on the assumption that friction is 


entirely negligible. Thus friction has no first-order effect. 
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The attenuation of the vertically-fronted wave as it travels up the pipe- 
line may be followed by putting ¢ = x in (3.11), which yields 

P Ue-tke, (3.25) 

Hence, on arrival for the first time at the upper end of the pipe-line, the 


wave has been reduced in magnitude to e ** times its initial value. 


4. Exit velocity falling uniformly to zero in time /, 

We now consider the consequences of shutting the valve in such a way 
that the velocity in the pipe-line close to the valve falls to zero at a uniform 
rate in time fy. Since the initial conditions are unchanged, (3.6) holds good; 
but of the two boundary conditions (3.7) and (3.8), the former needs 


alteration. For t > 0, x = 0, we take 
: [ ? U . 
€ = U(j—- — ) — = (t—t))?H(t—ty), (4.1) 
Zt a 
which satisfies the required conditions 
i=0 = = 6 
g n t 
f) t ty u : — i) (4,2 
ot ty 
t bo; u 0 
i at 4 ] 
hence é if ; (1 —e-Pte) _ | (4.3) 
\p' ty py 
From (3.6), (3.8), and (4.3) we then have 
z y a — P< | ge . 
c dé l V4. pte) Sinh q(] a) 8 kU+kUa (4.4 
dx t, p* cosh q p 
The substitution (3.10) will again be used; therefore we want the inverse 
transforms of anal 
u e- and l ani OF, 
p* p® 


The former can be obtained by means of the integration theorem from the 
inverse transform of gp-2e-%, and on integrating, between the limits ¢ and 
x, that part of (3.19) which lies within square brackets, we arrive at the 


expression 


ker 
ft - 
5) 13 Jey 19 J.2..2 Boyd _| [g(15 7 hy 1 122 
k (—8P--Bhka— jpk?a® — gy hea3 + kt gkat py hPa?) } 4 
_" 
f 9 [22 1 ].4>4 (33 3_[.2 2 
Et Soh v2 sigk*at4 kt(B— gh x*)-+ 


1 e242(19 


sigk*a*) | a | sigk*t4} H(t—z2). (4.5 
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If we apply the process to (3.20) we find 


at 7 —$ka—}k?a?+ kt(3+-4ka)}+ 


5) 


{7 beat eet yee] Hi x). (4.6) 


n both (4.5) and (4.6) the constants inside the curly brackets are equal 
n magnitude and opposite in sign. To avoid the difficulty which arises in 
inding a small difference between two large quantities, it is essential for 
yurposes Of computation to expand the exponential terms. Equations 
4.5) and (4.6) then become 





(4,2 


nverst 





1 | , ee 
Ira 3 he2x? 1h 3 eo raqasoK*x*) 
kt(1—ka+ py hk? Kart) + k2¢2(1 — ek? a? + gdsghketart) + 
Ke ( ‘8 Ka") 4 k*t"(sh4 sone") 
sink T slaok*t| H(t x), (4.7) 
Pats , or 
ka+3k22?— 1303+ 5 kat) + kt(1—ke+ £2’) 4 
| k 242/ | dk?) Ak | gagh tt | H(t x). (4.8) 


\s many terms as possible are given, subject to the proviso that none 
requires alteration if additional terms in the expansion of the Bessel func- 
tion are employed. Both (4.7) and (4.8) reduce to zero when ¢ is put equal 

To determine the other inverse transform we make use of the fact that, 


{ 
1 


j(p) is the transform of y(t)H(t—2), eP’%(p) is the transform of 
y(t—t,)H(t—t,—2). 


Hence in (4.7) or (4.8) (f—t¢,) is substituted for t, both in the unit function 
and in the remainder of the expression. In effect, no additional calculation 
of the terms typified by (4.7) is required; for a fixed position on the pipe- 
line, the values already found are used again with the times increased by fp. 
Thus the determination of P follows the same lines as in section 3, it being 
necessary to evaluate a succession of terms associated with H(t—z), 
H(t—t,—x), H(t—(2—2x)), H(t—t,—(2—z2)), ete. 
For a 0, the coefficient of H(t—2) in (4.7) and (4.8) is 


t+l1ki?—2h?8+-.... (4.9) 
Z 
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338 A. 
Thus as far as t = 2 or t = ty, whichever is the smaller, the initial pressure- 
rise does not follow a line that is exactly straight. Its gradient increases 
unless kt is very great, and then the approximation, which has been used, 
is known to be defective. 


5. Numerical results 

As examples of the methods set out in the two preceding sections, four 
tests have been considered which are fully described in the report (1) 
previously mentioned. For the conditions obtaining in these tests the 
pressure fluctuations at the valve (x 0) have been calculated. The 
details of the two pipe-lines and of the oil that they conveyed are given 
in Table 1. 


TABLE 1. Details of the pipe-lines and of the oils 








McCamey 4-in. | Wink—Monahans 6-in. 
Tests 23, 228 Tests 51, 50 
Pipe internal radius, a ft. 0-1677 0-2526 
Pipe length, / ft. : - ‘ 2232 108800 
Kinematic viscosity, v ft.2/sec. . 5-62 x 10-4 0-47 x 10-4 
Density, po lb./ft.% ; , : 56-4 52-2 
Wave velocity, c ft./sec. , : 4010 4000 
2.8) and (2.6) . 0-0892 0-160 


k, obtained from ( 


In the two experiments on the McCamey pipe-line the Reynolds number 
was below the critical, hence k as determined from (2.8) and (2.6) was 
employed. For the other pipe-line, however, the flow was turbulent and 
the value of k required adjustment. 


(i) Instantaneous closure 

For the McCamey pipe-line the value of k in Table 1 was checked by 
means of the results of test 23. From the well-known formula for the 
pressure drop in laminar flow it was deduced that at a steady velocity 
1-03 ft./sec. the drop should be 4-45 lb./sq. in. This is in agreement with 
the observation that the difference between the static and running pres- 
sures at the valve was 121—117 = 4 lb./sq. in. Accordingly, with the 
above value of k, the terms typified by (3.19) were evaluated at time 
intervals of 4, giving the ratio of the pressure P, (measured above the 
static pressure) to U. The calculations were pursued as far as ¢ = 6° 
where }kt was only 0-290, hence the approximation made in (3.18) should 
have no appreciable consequences. The results are shown by the broken 
line in Fig. 2, to which the ordinary rectangular diagram for frictionless 
motion has been added as a full line for purposes of comparison. After 
the initial jump equal to unity at f = 0, the broken line rises as more and 
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more of the oil is brought to rest and the friction head is reduced. Owing 
to the dissipation of energy which continues during the process of com- 
pression, the line fails to reach the ordinate 1-0 at t = 2 by an amount 


given by (3.24), which is here too small to be distinguishable in the figure. 
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Fic. 2. Pressure fluctuations at valve. Instantaneous closure. 
} 0 (frictionless) ; k 0-0892: . k 1-14. 


At t 


pressure is taken to be sufficient to prevent bouncing, the pressure suddenly 


2 the wave reflected from the open end arrives, and since the static 
becomes negative. At this instant the whole of the oil column is moving 
towards the reservoir; therefore, owing to friction, the pressure in it is 
As the 


tension wave retreats from the valve and brings the oil to rest, this friction 


sreater, i.e. the tension is less, at the valve than at the reservoir. 


effect is progressively reduced. Thus between ¢t = 2 and t = 4 the tension 
at the valve increases. The fluctuations continue in the same manner with 
slowly diminishing discontinuities. 

In contrast, the friction loss in the Wink—Monahans pipe-line was very 
great. In test 51 with a steady velocity 1-09 ft./sec. the loss was 56 Ib./sq. 
in., Whereas that given by the laminar flow formula is 7-87 lb./sq. in. The 


] 


value of k in Table | was therefore increased in the ratio of these figures, 


giving k |-14 as the value to be employed in the computations. The 
results are shown by the dash-dot line in Fig. 2. Here the friction is so 
large that the initial rise of unity is insufficient to bring the pressure up 
to the static pressure, and the maximum pressure, which occurs at t = 2, 
is 25 per cent. below the frictionless value. The ensuing sudden fall is too 
feeble to reduce the pressure below zero until a further lapse of time. The 
1-71, which is as far 


calculations were terminated at f 3 when skt 


as the approximation should be pushed. 
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(ii) Closure in time t, | 

As examples of the theory explained in section 4, the terms typified by 
(4.7) have been evaluated for both pipe-lines for ¢, 1. The value of x 
being fixed, a useful check on the computations occurs in the following way. 
When we are considering the term involving say H(t—(2—2)), (2—2) is 
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J 
Fic. 3. Pressure fluctuations at valve. Closure in time fy - 
k 0 (frictionless) : k = 00-0892; —- k 1-10. 


substituted for x inside the square brackets in (4.7) or (4.8). These ex- 
pressions vanish when ¢ is put equal to (2—2), as will be seen both from 
the algebra and from the obvious physical fact that the pressure change 
at the nose of the reflected wave is zero. 

For the McCamey pipe-line k = 0-0892 has again been taken. The result 
is shown in Fig. 3, whichis drawn with the same notation of lines as in Fig. 2. 


The initial pressure-rise steepens very slightly as it proceeds, and after 


¢ = 1 the pressure continues to increase as in Fig. 2. The effect of friction 
upon the maximum pressure attained at tf = 2 is now more pronounced, 


and at this point a visible difference exists between the full and the broken 
lines. 

For the Wink—Monahans pipe-line it was desired to make a comparison 
with test 50. In this the initial steady velocity was 1-05 ft./sec., which was 
somewhat less than in test 51. The corresponding value of k, obtained by 
the same method as before, was 1-10, and this was used in the calculations, 
yielding the dash-dot line in Fig. 3. Although k is less, friction has again 
a larger effect, amounting to a reduction of 41 per cent. at ¢ = 2. 

(iii) Comparison with the graphical method 

In the report the results of the graphical method are plotted as isolated 


points on reproductions of the pressure-gauge readings. For tests 23 and 
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22B on the McCamey pipe-line and tests 51 and 50 on the Wink—Monahans 
pipe-line these points have been transferred to Figs. 4, 5, 6, and 7, where 
they are shown surrounded by circles. To these figures the results (with 
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Fig. 4. McCamey pipe-line (test 23). Instantaneous closure at 1-03 ft./sec. 
2 < by graphical method ; by calculation. 
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test 228). Closure in time l/c at 1-03 ft./sec. 


Fic. 5. McCamey pip 
. *) by gray hical method ; by calculation. 


| triction) displayed in Figs. 2 and 3 have been added in full lines after 
transformation to the dimensional quantities appropriate to the tests. In 
Fig. 4 it appears that in the graphical method no attempt was made to 
obtain the sloping part of the diagram; moreover, one end of each dis 
continuity is omitted. Nevertheless, the two methods agree in yielding 


the same maximum pressure, and this is true also of Fig. 5. 
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For the Wink—Monahans pipe-line (Figs. 6 and 7) with its heavy friction, 
the discrepancies are more serious, and again the discontinuities resulting 
from sudden closure are blurred in the graphical method. An explanation 
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Fic. 6. Wink—Monahans pipe-line (test 51). Instantaneous closure 
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Fic. 7. Wink—Monahans pipe-line (test 50). 
at 1-05 ft./sec. 


>) by graphical method ; 


Closure in time l/c 
by calculation. 





‘an be put forward of the fact that the graphical method yields a larger 
maximum pressure, the difference being more marked when closure is 
gradual. The friction—loss relation which was employed does not seem t¢ 
be explicitly stated in the report, but to judge from a worked example it 
almost certainly was taken to vary as the square of the velocity throughout 
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every test. On the other hand, a loss proportional simply to the velocity 
is assumed in this paper. Thus, as the velocity fails after the closure of the 
valve, the friction loss in the graphical method diminishes more rapidly, 
leading to a higher maximum pressure. Clearly this effect should be more 
marked when the oil is brought gradually to rest; no doubt it is too slight 
to be noticeable in the McCamey tests, in which friction was so small. It 
would be out of place in a computational paper to discuss the experimental 
observations in detail. But it may be remarked that, whereas in tests 2: 


and 228 the graphical method satisfactorily predicted the maximum 


pressure, in tests 51 and 50 it led to a considerable over-estimate and the 
other method is more accurate. 

The report also gives details of a test on the Wink—Monahans pipe-line, 
in which the velocity was raised to 4-38 ft./sec. by means of a pump. Under 
these severe conditions the effective value of k is so large that the approxi- 
mation (3.18) cannot provide an accurate result for the maximum pressure. 
No doubt this case could be met by employing more terms, but only at 
the cost of considerable labour. Evidently it will not be possible to claim 
that the general problem with laminar flow has been solved completely 
until tables of the integral in (3.11) are available. 
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The problem treated in this paper has been examined also by M. Lupwie (Amer. 
Soc. Mech. Eng., Preprint No. 50-SA—26, 1950). 











ON THE HALF-PLANE DIFFRACTION PROBLEM 


By F. G. FRIEDLANDER (The University, Manchester) 
[Received 22 November 1950] 


SUMMARY 
The diffraction of an arbitrary two-dimensional disturbance by a semi-infinite 
plane screen is considered. It is shown that the solution, at all points of space, can 
be obtained by a modification of the method developed by Hadamard for Cauchy's 
problem. The diffraction of the disturbance represented by Hadamard’s ‘elementary 
solution’ is discussed, and a remarkably simple expression for the diffracted distur. 
bance is found; its Laplace transform is shown to agree with the Green’s functions 
derived from Sommerfeld’s two-valued solutions of the wave equation. 
1. THE diffraction of two-dimensional waves—both infinite harmonic wave 
trains and aperiodic disturbances—by a semi-infinite plane screen has been 
treated by various authors since Sommerfeld derived an exact solution for 
the case of an incident plane harmonic wave train. Various methods of 
attack, ranging from Sommerfeld’s many-valued solutions of the wave 
equation to the reduction of the problem to an integral equation of the 
Wiener—Hopf type, have been employed. In the present paper it is shown 
that the problem can be solved by a modification of the method developed 
by Hadamard for Cauchy’s problem. 
We shall consider the two-dimensional wave equation 
Cu =O tu 
—— ——-__. —= f(x, y, t), a 
a= ox* day* 
where x, y are Cartesian coordinates, ¢ is the time (on a suitable scale), 
f(x, y, t) is a known function, and wu is a variable whose physical significance 
depends on the problem under consideration; e.g. for sound waves u maj 
be the pressure, and for electromagnetic waves it may be a component of 
the Hertz vector. It is more convenient to consider two auxiliary prob 
lems, which may be formulated as follows. 
Let p(x,t) and q(x,t) be defined for « > 0, and be identically zero fo1 
t < o(x), where |o'(2)) < 1;+ they must also satisfy certain conditions 01 
regularity which will be mentioned later. In the ‘space-time’ (x, y,t) the 
union of the influence domains of the points (a, 0, o(a)) is bounded by a 
characteristic C, which is the envelope of the characteristic half-cones 
t = o(x’)+{(x—2')PP+y*}* (O< a’ < om). 
p(x, t) and g(a, t) may be considered, jointly, as a disturbance propagating along th 


half-plane. The condition |o’(«) 1 means that the velocity of the front of this disturbance 
does not exceed the velocity of wave propagation. 


[Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 3 (1951)) 
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It is easily shown that if a point (x,y,t) is abovet C,, then the points 
y (x’, y’, t') which lie above C,, and for which 
! < t—{(x—2’ P+ (y—y'}*}, 

are a bounded set of points. We then consider 
PRoBLEM I. 1'o find a function u,(x, y, t), regular except on y = 0, x > 0, 
satisfying (1.1) and 








ni . ¢ Uy . > 
(a) lim q(x,t) (x > 0); (1.2) 
0 CY 
a ‘ 
ental (b) Uy 0 on C\: (1.3) 
dis 
| 
ict : ou ou 9 9 
(c) t hig 1.0 as a4 y* —> 0. (1.4) 
Ox oy 
Wave = . ’ 
ProsBLEM II. 70 find a function u,(x, y, t), regular except on y = 0, x > 0, 
Ss DePT | . - . “5 . . 
| satisfying | ) and 
on for | S@# fying (1.1) ane 
rds of (a) lim U5 pl(x,t) (ax 0); (1.5) 
y=0 > 
wave 
of the (b) condition (1.3); (c) condition (1.4). 
show By a trivial modification of the proof of the well-known uniqueness 
sloped theorem, it can be shown that the solutions of both these problems, if they 
exist, are unique (1). (The condition (1.4) excludes certain diverging 
evlindrical waves, infinite at x = y = 0.) 
Now let u(x, y,t) be a solution of (1.1) which vanishes on and below a 
characteristic C, which intersects (y = 0, x > 0) in a curve t = o(x) (then 
necessarily |o’(x) |). If we define C, by means of this curve as above, 
ee and put 
icance ou » 
| Pv t) | Uo] 0 q(x, t) | | (x Q), (1.6) 
“4 ma ly 
' cy y=0 
ent , 
then the functions 
pro! u U, Uy. U4 Uo Us (1.7) 
ro f satisfy the boundary conditions 
ons l¢ Us | , 
aa i— 0 iMlug = 0 (2 > 0). (1.8) 
, 1 cy 1=0 
lb ™ 
- hus they are the solutions of the diffraction problems for the semi-infinite 
screen (y = 0, a 0), the incident wave being uw», and the boundary con- 
NtIONS CU/CY 0 and u 0 on the screen, respectively. 
u I 1X thought of as drawn vertically upwards, so that ‘above’ is equivalent 


is equivalent to ‘earlier’. 
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2. Let S be a closed surface of space-time, (/,,/,,/)) the direction cosines ap} 
of its inward normal, and du/dv the transversal derivative defined by we 
du ou ou ou 
=1,——1,—— 1, —. (2.1) 
dv ot Oa oy 
[f wu satisfies (1.1) and v satisfies the homogeneous wave equation 
on Czy On 
—. —-—___.__. = Q), (2,2) 
ol =6ox* = 6oy* 
then, by a modification of the divergence theorem, 
i dv du\ 5. ? 
| (u —— 9 —1 05 = = [ fe dT’, (2.3) 
J \ dv dv, 5 
where dS is the surface element of S, dT’ the volume element dxdydt, the 
surface integral is taken over S, and the volume integral over the interior 
of S. This equation is in the first instance only valid if u and v are suffi 
ciently regular in and on S; but Hadamard has shown that it can be used 
even when S contains portions on which, for example, v is ‘fractionally 
infinite, provided that the integrals are interpreted as ‘finite parts’ of 
divergent integrals, in the sense defined by Hadamard (2). 
Let (x, y, t) be a fixed point of space-time situated above C,, and (2’, y’,t’ 
be a variable point. The elementary solution of (2.2) is 
'\2 12 1)2\—-4 ‘ 
v {(t—t’)?—(4—2'’)?—(y—y’)*}-#. (2.4 
It can be considered as a ‘retrograde’ wave whose wave fronts constitute 
the characteristic C, 
F aut S(> »/\2 "\2U4 95 
t= t—{(x—2')°+(y—y')*f'. (2.5 
With it, a retrograde diffracted wave is associated whose wave fronts 
constitute a characteristic C,, 
t’ = t—(x*?+y?)*§—(x’2+y’2)}, (2.6 
(see Fig. 1). Put 
x = rcos8, y = rsin®@, x = cost’, y =r sind’, (2.7 res 
the 
a ff’ “\2_] , 2.4 “ -2_|_ »’2 ory’ ang , 14 9 
R = {(v’—x)?+(y’—y)?}! = {r?+7r2—2rr’ cos(0’—@)}}. (2.8 des 
We shall take, in the first instance, y > 0, so that 0 < 6 < 7. The equa- 
tions of C and of C, become Th 
' the 
t’ = t—R, t’ = t—r'—+ (2.9 
pre 
respectively, and it is obvious that, for 0 < @ < 7, C, is below C except 
for 6’ = 0+-7, when they touch along a generator. 
la Bl y , ° , , s ‘) Th 
3. The surfaces C, C,, for 0 < 0’ < 64 nm, together with (y 0,2 ov 


and ©,, delimit a bounded domain of space-time to which we shall nov 
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apply (2.3), taking v to be the elementary solution (2.4). In order to do so, 
we first introduce two auxiliary characteristics C’, Cj, given by 


t’ t—(R?+-e)}, t’ = t—(r?+-€)?—r’ (3.1) 


At 
























































Fic. 1. 


respectively, where ¢ is afterwards made to tend to zero, and apply (2.3) to 


the domain bounded by ©’, C’, (y 0, x > 0), C,, and a surface S 
defined by 


” 
;, A ‘ 
t i— », R< (7° -¢)? (y > eé*). 
The contribution of this surface S, to the surface integral as « > 0 first and 
then » + 0 is evaluated exactly as in Hadamard’s treatment of Cauchy’s 
problem and turns out to be (neglecting fractional infinities in e€) 

27ru(2, y, t). (3.2) 
Che surface C, makes no contribution to the surface integral; for in either 
problem, « = 0 on ©,, and as C, is a characteristic, the transversal to O, 
is also tangential to C,, so that also du/dv = 0 on C,. 
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On (y = 0, x > 0), the transversal derivative is —0/éy and so its con- | gybs 
tribution to the surface integral is, as « + 0, 


* 


> 


“i 


oY 


ov ou 
+2 
cy oy 
where the asterisk indicates that this is to be interpreted as a ‘finite part’, 
and the domain of integration is bounded by 





| dxdt’, and | 
y’=0+0 


t’ = o(2’), t’ = t—r—2’, t’ = t—{(x—2’)?+- yy, wher 
Puttin li (a, t) lj ou (at) ( 0) (33 The 
u g im u (a, t), im — = q(x, x > 0), a3 
i y= +0 f y= +0 0Y , tend 
and R, = {(a—2’)?+ y?}, (3.4) 


this can be written as 


dt’, (3.5) 


ae ee 3.5) |" 
J UG’ P— Ro} (tt)? — Bo}}) 





where the upper limit oo has been used for convenience; as in fact p(2’, 1’) 
and q(x’,t’) vanish for t’ < o(x’), the domain of integration is finite. 


“a ee canes cae sal TI 
rhe contribution of C’ is of the form A(e)e-?, where A(e) tends to a limit tal 
ye : esta 

as e > 0, and hence its finite part is zero. 
into 
tend 


4. It only remains to consider the contribution from C),. On C’, we can 


dS ad =— —?7’ was _ “) 
dv 





write 
et’ er} 

Also, if for the moment we distinguish functions evaluated on C’, by en- 

closing them in square brackets, it follows from (3.1) that 


C ou ou C ov Ov 
alt] a E ' ap lt in l= (FI. wher 





’ , . " . ° : U 
Hence the contribution from C, to the surface integral is | 
Qin r whe 
7. ‘ ; a 
r= | ao’ | \(uj Spo}—[o] Sql ar’, (4.1 
I *o er 
0 0 
where the upper limit oo has again been introduced for convenience, | 
vanishing when 7’ is so large that the corresponding point of C’, lies below 
C,. In view of this, and of condition (1.4), a partial integration transforms 
(4.1) into P m 
O+n Thu 


r= [a | [ujfzr fo) +(e} dr’, 
. P cr 
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Substituting the value of |v], which by (2.4) and (3.1) is 


[v] fet 2r’(r?+-€)*+ 2rr’ cos(d—8@’)! 


ind putting @ 7+6—d, this becomes 
Le 
” eti(v’, 7+0—d¢d, t—(v?+€)?—?’ ? 
; | dd | eee (4.2) 
} fet 2r’(r?2+-€)?—2rr’ cos dh}! 
0 0 
where u(r, 6, t) u(r cos 0, r sin 8, t). (4.3) 
[The behaviour of the kernel of the integral in (4.2) is singular as e > 0; it 
tends to zero with e if both -Qand d ~ 0, but at the same time 
a t-A 
‘ F € dd r } 1 
lim : ——— — = 1° = (4.4) 
e-o0 J {e+ 27'(r?+€)?—2rr’ cosd}? r r+r 
0 


Hence it can be shown—assuming that uw is continuous—that 


* alr’, 7+, —r’)[r\4 , 
Me. (| Sl (2) ar. (4.5) 
e=0 T—+—-f" r 


The detailed proof—which is not difficult—will be omitted, but we shall 
establish (4.4). If the range of integration in the integral in (4.4) is divided 


into (0, n) and (y,7+-6@), where » 0, then the contribution from (7, 7-+-@) 
tends to zero with « for fixed ». We need therefore only consider 
f 
edd 
J | 
J 4 2) -¢€) 2rr’ cos hb} 
t] 
[ € dias 
: (1 tes")? t« 2r' (r?2 +e)? —2rr’- rr’ ys? 
0 
where cosd Lj? 
Clearly J J J* sec 3, 
where 
dis 
J* mse 
P 2r (r=+-« 27) Yr) ‘yb?} 
€ 2 sin by 
e+27'(r?+e) 2 fet 2r'(r2+-e)t—2rr’ 477 sin? 3y}4 
= ' ] r\4 
: Thus lim J* — ie 
é ) r r r | 


and letting finally - 0, we obtain (4.4). 
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Collecting our results, (3.2), (3.5), and (4.5), we obtain now 


- fia (q(x', t’) , yp(xt)| it’ — 


da’ | x 


277ru(x, y, t) : ft 
J | Yo Yo J 





ere a ne ee. rn tak wet ae ae 
| u(r", @ 7, t —- \e) dr’ 4 | | | file y st) dx'dy'dt’, (4.6) 
! RY r 


r 


where y = (t—t’)?— R?, Yo = (t—t’)?— RR, (4.7) 


ir er<cek 660 < 22, 


and f, is equal to f above C,, and zero below C\. 


5. Again, let (2.3) be applied to a domain bounded by ©,, C, and by 
6’ = n, 0’ = 2x—7, v being still the elementary solution (2.4). As » +0, 
this domain tends to Y,. The contributions of 6’ = », 0’ = 27—7 tend, 
as 7 > 0, to terms of'the same type as (3.5), with p, q replaced, respec- 
tively, by 


ou|y= +9 = 
| ‘ (5.1 


p*(z,t) = [uf-*.. q* (x, t) E 
CY | 


y=—0 


and the limits of the inner integral replaced by (—oo,t—r—za’). As (x,y, 
is outside the domain, there is no term corresponding to (3.2). As for the 
contribution of C,, it is calculated as above: but as the sense of the 
transversal must be reversed, its sign must be changed, and as the integral 
corresponding to (4.2) must be taken over (0, 27), the result (4.5) must be 
multiplied by 2. Hence after division by 2 the resulting equation is 


v 


| L 
aye 
. ‘ 70 70 
f 


omy fan f (Met, wret ay. 


y3 


, u(r’, 0-2, t—r—r’)[r\4 k tt hie a € rr ryt =D 
| ur is = eee (5) dr’ 4 | | | file oy ai dy dt’ (5.2 


a r r 


where the volume integral is taken over t’ < t—r—7?’, 0 < 6’ < 27. 

If this equation is added to (4.6), the contributions from C,, cancel, and 
we obtain the following equation, which must necessarily be satisfied by 
any solution of the wave equation which is sufficiently regular to ensure 
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the validity of the various limiting operations by which it has been 
obtained: 


7. 7 


> y’ ) } a 4 ‘¢ a , : , , 
2ru(a, ¥, t) hit a4 0) dx'dy'dt' ! } | | | filx’ y é ) dx dy dt 
ova x ag aaa i 


0 d 
% * 
rr fy 4 7! r(x’. , “pn * x’. , * x’. , : , 
~~) es eer | | (a8), POON awat 
l ¥% YO J. Yo Y0 
/ g; 
(5.3) 
where Y,, GY, F SF, are defined by 
D, t—r—) t t—R, 0 ’ a+6; ) 
= t t 0 27; 
d X r—? (5.4) 


This result holds for y 0; if y < 0, the sign of the surface integrals must 
be reversed, and in the definition of %, 6 and 6’ must be replaced by 
2r—O0, 2rn—O respectively. 


Now let D y-* {(t—t')?—(xa—2a’)?—(y+y’)?t-3 (5.5) 


i 
/ 


be the elementary solution relative to the point (7, —y, t) which is the image 
of (a, y, t) in y 0, put 


R ((~—a’)*+-(y+y’)?} fy21 »’2 rr’ cos(O+-0’) 4, (5.6) 
and let -. be defined by 
G t—r—) <t—R, 0< 06 < a—8, (5.7) 


lor y > 0 (see Fig. 2); if y < 0, 6 and @’ must again be replaced by 27—8, 
27—@' in the definition of ZY. Applying (2.3) to ZY, with 6 instead of »v, 
md to Y, with 36 instead of v, and adding the resulting equations, we 
obtain, after a calculation similar to the one which has just been carried 
out, that 


2 rn Pp 


A | ae . ad. da'dy'dt 2 | ] [ Ay: ) dx'dy'dt’ 
ar ie Jvc r 


ay? 


70 70 


+ * 

. . a yr’ t’) ) } ind a*(a’. t’ ¥*(a'. , 

He Wee ED antae—y | | PRESS) _ aE arag. 
vs Yo J « | Yo 


(5.8) 


the left-hand side is zero as neither Y, nor Y, contains the point (x, —y, t). 
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6. We can now deduce the solutions of problems I and IT from (5.3) and 
(5.8). In problem I, q is given and so, by (5.1), g* = 0. Hence, adding (5.3 


(5.3) and (5.8), we obtain 
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where V, is defined by Lip 


.- ey a ae eR SS re (x, 
L=—= mG; yi=-—+=— mG; 4,= 3(- +=} m DY, ons 
y? y 7 . 2\yi a 

(6.2 A 


7 . . . . . . > : 31] exXa 
J, being the domain of points belonging to ZY, but not to Z,, and it will xa 
- 0 below C,. — 


be recalled that f, = f above C, and f, 
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In problem II, p is given and hence p* 0; subtracting (5.8) from 
5.3). we obtain therefore 


J “rc (x’. t’) ar 
U(X. Y, t : | f da'dt' 4 
, V(x, y.tix’ yt) f(a’. y',t') da'dy'dt’, (6.3) 
where 
| l l : 1/1] ] 
y, in & A —~ nZY Vp ( _— | in Y, 
7 7 “ 


(6.4) 


The solutions of problems I and II, if they exist, are necessarily of the 
forms (6.1), (6.3) respectively. To make the argument rigorous, we should 
now establish the existence of these solutions by verifying that the 
expressions (6.1), (6.3) for uw, and uy satisfy the wave equation and the 
respective boundary conditions. The calculations involved are lengthy, 
and of some complexity, and will be omitted. If the homogeneous wave 
equation is considered, so that f = 0 and the volume integrals in (6.1), 
6.3) are absent, no particular difficulty is encountered. To verify that the 


wave equation is, for example, satisfied by 


=— | | = - ) da’ dt’ 
wT de Yo 
Ai 

it is permissible, according to Hadamard, to differentiate under the integral 
sign, provided that the resulting integrals are interpreted as ‘finite parts’, 
and that any fractional infinities arising from differentiation with respect 
to the variable limits are omitted. (The expressions for the second-order 
derivatives also contain finite line integrals; these, however, cancel on 
substitution in the wave equation.) This process is valid only if the finite 
parts of divergent integrals occurring in the formal expressions for the 
second-order derivatives have a meaning. A sufficient condition for this 
is that g should possess partial derivatives of the first order satisfying 
Lipschitz conditions. But if g is, for example, singular on a curve in the 
(x,t)-plane, the expression for u, may still be valid; such cases are best 
examined as they arise 

\s for the boundary condition at y = 0, x > 0, the verification is 
exactly as in the case of the boundary y = 0, o<e< oo. It is 


analogous to the treatment of a double layer in the theory of potential. 


5092.1 


Aa 
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The volume integral is more troublesome. The singularity near 2’ = g, 
y’ = y, t’ = t occurs in the solution of Cauchy’s problem, and is dealt 


with in Hadamard’s book. However, formal differentiation gives integrals 
with singularities at 6’ = 70, t' = t—r—r’,r’ > 0, to which the method 
of the ‘finite part’ cannot be applied directly. I have not carried out the 
calculation in this case, but there is little doubt that it would lead to the 
requisite verification. 

The solution (6.3) can be considered in a similar manner. It must be 
noted, however, that as the surface integral is now a ‘finite part’, a suffi- 
cient condition for the validity of the solution is that the partial derivatives 
of the second order of p should exist and satisfy Lipschitz conditions. 


7. A little consideration will show that the process by which (6.1) and 
(6.3) have been obtained is equivalent to taking the function v in (2.3) to 
be V,, defined by (6.2), or K, defined by (6.4), respectively. These functions 
therefore play the role of Green’s functions for these problems, and by 
reciprocity the solutions of the associated diffraction problems for an 
incident wave generated by an instantaneous point source can be obtained. 
In fact, suppose that p = q ), and let f(x, y,t) tend to the singular 
function 
f(x,y, t) = 8(a—2_)d(y—y)8(8), (7.1) 
where 6 denotes Dirac’s 5-function. The limiting values of (6.1) and 
(6.3) are then the solutions of the diffraction problems with an incident 


wave 
Ug {t? (x a)" (Y -Yo)?} t+ (> ((%—2X9)” +(Y—Yo)"s") ) | (7.2) 
0 (t < {(@ aX)” +(Y¥—Yo) )?}*) 
for the boundary conditions éu/éy = 0 and u 0,aty = 0, x > YO, respec 


tively. Let, for example, y, > 0 and put 


p {(a—ay)?+(y—Yo)*}}, p {(x—a)?+(y+Yo)?}? ) 
19 COS Ay, Yo rosin, (0 <6 <7) 


Xo 


(7.4) 


Let the domains 7°, 7°, 7 be defined by 
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Then it follows from (6.1) and (6.3) that the solutions of these ‘elementary’ 
diffraction problems are 


u u, in Z u UjtU, in D; | 
0 oO r (7 6) 
u 3(Up+%,) in F; u = 0 at all other points, J 


with the upper sign for problem I and the lower sign for problem II. (If 
the origin in Figs. 1 and 2 is moved to (0,0,t), the sense of the ¢-axis is 
reversed, and 2, y are replaced by 2», Yo, then these figures represent 
0 fon 0 C70 


P+ and Zi, and GZ, G9, respectively.) 
These results are remarkable, both because of their simplicity, and 


because of the discontinuity on the ‘diffraction characteristic’ 


t Yor? 


In either case the discontinuity is compatible with the wave equation, or 
rather with the generalized continuity equation 


Onn 
ds 0. 


(Lv 


. 


which follows from (2.1) by taking v = 1, f = 0, as the square of the saltus 
satisfies the intensity law of geometrical optics (3). Mathematically, the 
discontinuity is probably due to the singular behaviour of the elementary 
solution (6.5) on its wave front. Its physical meaning is obscure, but it 
must be remembered that (6.5) represents an instantaneous infinitely long 
line source. For a source which is not instantaneous, the discontinuity 


disappear S 


8. Finally, it may be of some interest to compare the results obtained 
here with some of the results obtained by different methods. In the first 
place, consider (7.6). It is to be expected that the Laplace transform of 


this solution is the Green's function of the half-plane for the equation 


eu al a 
: = sil , (8.1) 
Ox* cy* 
Let (x,y) be a point ol the shadow, so that 7-4 G, <0 < 2a. Then, by 
7 5 (7 8) 
"u (Up Uo f } To); u 0 (f r+ fo). (8.2) 


this function is 
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Putting /? = p?+-* in the first, and # = p?+v? in the second of these 
integrals, we obtain the expression 
x x 


; l . s(p? ip2yt l . (57 +2) 
U =- cay“? = dv, (8.3) 
Z N (p*--P Z J yee 
2(rro)|cos{(8—4)/2 2(rro)\cos{(8 + A)/2} 

with the upper sign again for the boundary condition éU /éy = 0 ony = 0, 
x > 0, and the lower sign for the boundary condition U = 0 on y = 0, 
x > 0. This, however, is equivalent to the Green’s functions deduced from 
Sommerfeld’s integral (4, equation (42)); it is easily verified that at points 
not in the shadow, the corresponding expressions are obtained by taking 
the Laplace transform of (7.6). 


Next, let us consider the homogeneous wave equation, i.e. take f = 0. 





Then, if we let y > +0 in (6.1) we obtain, for x > 0, 
f t r sf 
: ey : r q(x’, t’) dt’ 
lim u, = 3 | dx . a = ‘= — (8.4) 
yin tt 7. J) {tt P= (@— 23} 
0 t—x-—a2 


which is equivalent to a result which has been obtained by Evvard (5) and 
Ward (6). 

In the homogeneous case, u,(x, —y,t) satisfies the same boundary 
conditions as u,(2,y,t), as well as the wave equation, and so by the 


uniqueness theorem 
1 U,(x, y, t) U,(x, —y, t), 


so that uw, is an odd function of y. In particular, it vanishes on y = 0, 
x < 0, and so by (1.7) the function u,, which is the solution of the diffrac- 
tion problem for an incident wave uw, with the boundary condition 
éu/ey = 0 on the screen, satisfies the boundary corditions 


u u (y O, & <= 0); Cujdy = 0 (y= 0, x > 0). (8.5) 


0 
Similarly, w, is an even function of y and its derivative vanishes on 
y = 0,x < 0. Hence, if the wave-front characteristic of wu, is of a suitable 
type, wu, may be defined, for example, in y < 0 as the solution of a problem 
of type I, with w given, and equal to uv», on y = 0, x < 0. This formula- 
tion of the diffraction problem has been used by Fox (7), who has obtained 
the value of uw on y 0, « > 0 from an integral equation. His result can 
also be deduced from (6.3), by considering the limiting value of the surface 
integral as y > 0 for x < 0 (since (6.3) applies to problem II in the un- 
modified form). 

Let us first calculate for this purpose 
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By definition, this is the limit of 
y : dr A 
| (7? aT 
Ro +« 
as e-> 0, where A is chosen so as to cancel the fractional infinity. From 
this it follows that 


y rx’ l/r—ax\h r+a 
J er <= if re 59 
Thy 4 2X (2 x)s W\ 22 J (X—z )*--9y" 


and hence, as y + 0 and so r > —2 (for x < 0), 
L/ja\\4 l i 
lim J | ‘ae : (8.7) 
0 T\ © = Zz 
Now we can write (6.3) for f 0 as 
L t—Ro 


, | p(a’,t—Ry)J da y | da’ | p(x ,t ~ Fo) pit ) dt’. 
{(t—t’)?— Re 


Assuming that p satisfies a Lipschitz condition, and is bounded, it is easily 
shown that the second integral tends to zero with y, and that the passage 
to the limit under the integral sign in the first term is legitimate, so that 


l f p(x’, t+a—2’)/\x\\3 , 
lim wv, mods a ( =) dx’ (« < 0), (8.8) 
MY 


0 7T Hy 


which is equivalent to the result obtained by Fox (7, Appendix B). 
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A NEW RELAXATION TREATMENT OF FLOW 
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SUMMARY 

In this paper methods are developed for the calculation of incompressible and 
compressible flow under conditions of axial symmetry. Relaxation is used to solve 
new forms of the flow equations. In these equations ¢, the velocity potential, and 
us, Stokes’s stream function, for the incompressible flow, are the independent variables. 

The radial distance, 7, is the dependent variable in the incompressible flow equa- 
tion, and this feature renders the method especially suitable for the design of axially 
symmetric ducts. A typical design problem is calculated as an example at the end 
of the paper. 

Log(1/q), where q is the velocity, is the dependent variable in the compressible flow 
equation, which can only be solved after the incompressible flow solution has been 
found. 


Nomenclature 
r radial distance measured from axis of symmetry. 
x distance measured parallel to axis of symmetry. 
u,v velocity components in the x and r directions respectively. 
(q,8) velocity vector in polar coordinates. 
(d,s) the incompressible flow plane, in which ¢ is the velocity potential, 
and % is Stokes’s stream function. 
n,m mesh size in ¢ and %& directions respectively. 
f ”. 
L log(1 q)- 
density. 
M local Mach number. 
local velocity of sound. 


> 
o 


~ 
~ 


¢ vorticity. 


be 


angle between the incompressible flow vector and the compressible 


flow vector. 


1. Introduction 

THis paper extends the methods of references 1 and 2 to the flow of an 
inviscid fluid with axial symmetry. Even for incompressible flow, the 
equation of motion is too difficult to solve except by a relaxation treat 
ment. Relaxation and ‘squaring’ treatments of similar problems of axial 
symmetry have been given (3, 4) but they have all been based on equations 


{Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 3 (1951)] 
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in which 2 and r have been the independent variables. This has the dis- 
advantage of requiring ‘irregular stars’ in the relaxation mesh near curved 
boundaries. The method of the present paper avoids this difficulty by 
taking # and &% as the independent variables. For the incompressible flow 
is taken to be the dependent variable, a procedure which permits the direct 
design of ducts 

The design problem is essentially the calculation of r on the duct bound- 
ary such that a specified pressure distribution is achieved. When working 
in the (vw, r)-plane with ¢ (or ys) as the dependent variable, the relaxation 
pattern involves r, which is. of course, unknown on the boundary in the 
design problem. (However, a ‘trial and error’ method of design in the 
x.r)-plane has been developed (5).) Interchanging the roles of (#,r) and 
é.u) avoids the difficulty. and the design can be carried out directly. As 
in example of the method a design problem is worked out in section 7. 

For compressible flow L (log(1/q)) is taken to be the dependent variable 
in the (d.y)-plane. The method is very similar to that developed previously 


for two-dimensional flow (2 


2. Incompressible flow 
By definition we have 
l 
db, u ys. d, v - yb .. (1) 


} r 


while Laplace's equation for ¢ and % can be written 


us d, 
Wf; i). d,., 4 E all Dp» 0. (2) 
} > 
where suffixes denote partial derivatives. Now 
f / 
CLD, Ww) 
ede / ee a 2 
us dus, r(u=+-v-) rq-. 
Ax.) 
ind therefore 
Ys, d, ub, ; db, 9 
»? . 9? 'dt 9° ty 9? (») 
rq- rq- rq : rq- 
where suffixes denote partial derivatives. 
From (1) and (3) we find 
"4 
r r? Cy PS (4) 
1 - 
: ( l or c cr ~ 
ind hence | r = @, (9) 
@ \) Ch ods os 


— ( l ox ( CXL 
and | r 0. (6) 
Cdh\r Cd c us ( ub 
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If we write ju? (7) 
—* éflef\. ef 

equation (5) becomes — “ + J 0. (8) 
cb\f ed Oys? 


Equation (6) is of no value as it involves both dependent variables, but 
equation (8) is relatively simple. A relaxation treatment of (8) is given 
in the next section. 
From (4) it follows that 
dx = rr, dd—(rg/r) df 
shy dd—(fa/f) dp}. (9) 
which enables x = 2(¢,y) to be calculated directly from the f = /(d.d) 
solution of equation (8). 
Further, from (1) and (3), 
q (75 + r?(ry)?)-*, ? = tan (4) 


ry 


i.e. q = 2AffF+(fy)*). 8 tan-1( Ie } (10 
wh 


} 
and so the velocity vector (q,@) can be directly calculated from the / field. 


3. Relaxation solution 
Fig. 1 shows four typical rectangles of the 4, 4 mesh. The points a. b. ¢, 
and d bisect the intervals 0-2, 0-3, 0—4, and 0-1 respectively. The differ- 
ence equation corresponding to (8) can be writtent 
(“) (2 fy_fe hs) 
(0 Sa fe 


For an arbitrary arrangement of the f field, X (residual of the relaxation 


+ (fotSa—2fo) = Xo = 9. (11) 


process) will not generally be zero. From (11) we find the ‘influence 


é Xo 7 c Xo lo. ("")( ' x . ) 

Chas Cfo | n}\fa fod) . (12) 
OX, [m\? /, OX, m\* | . j : 
of (7) [te fs (") |% 
of _of 


C cu 


coefficients’ 


For most practical problems and so it is sufficient to interpolate 


for f; and f, using 
fi 2(f3+So): fa = h+So)- 

Since f is usually large compared with cf/é¢ and ¢?f/éd¢?, a close approxi- 

ef 

xp? 


+ See Appendix for an alternative method. 


mation to equation (8) is 
0 
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c ‘ 
Le. u const., on @ const., (13) 
Ow 


and the corresponding difference equation is 


¥ 























Fic. 1. 


From the last of equations (3) we find that equation (13) implies that 
uq° = const., on d = const. Now the linear theory of axial symmetric 
flow assumes that v = 0,q¢ = u = const., on ¢ (i.e. in this case, x) = const. 
and so the approximation (13) is equivalent to assuming the linear theory, 
and is a useful starting assumption for the f field. 


4. Boundary conditions 

Normally either r or q will be specified on a boundary as a function of x. 
We shall assume for the moment that we have been able to find q or r as 
functions of ¢. Then with 7 r(d) there is no difficulty, while with 
q = 9(¢) we proceed as follows. 

Suppose the line 3-0-1 in Fig. 1 is the boundary, then from (10) 


6m* (=) Sa) (fo—fa)?: (14) 
I ii Jo va 


bs f, f, of4m* eeranie (15) 


| g@ n 


4fy 
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< 


Using this in (11) we have for a boundary point 
, _ (m\(fi—So , fs—So\ op . oft [m\Pi—Ss)"\3 oy 
Xo = |— ee fe lt ig ee a fo. 
n Ja to | la n | 
(16) 


and since usually (/,;—/s)?/4f) is 4m?qz, the approximate influence 
coefficients 
0X ) eXo 


of. Gh, 7 (7) [fe = ~ (r) [to 
e-em 


can be used. Any errors introduced by this approximation are detected 
when (16) is used to recalculate residuals as a check at the end of the 
calculation. 

On lines ¢ = const. at a sufficient distance upstream and downstream 
the linear theory will enable fairly accurate boundary conditions to be 
established. 

Now it remains to show how we can arrive at q q(h) from q = q(x 
Using the linear theory we can deduce an approximate ¢ = (2), and hence 
an approximate g = q(¢) to start with. The f field is then relaxed with 
these approximate boundary conditions. Next we integrate along the 
boundary, const., using a difference equation form of (9), i.e., referring 


to Fig. 1, 


—_—_ 
~ 
—— 
t= 
te 
he 
~~ 
— 


Substituting the value of f, from (11) we have 


—_ M\fop  ¢) Alm (= to, ts—So\| 9 
oa: a (“\lr Ja) (=) 7 Tf, }}: (I 


from which we find immediately x x(). 
Using this in g = g(a) we find a new and more accurate relation g = q(¢ 
the relaxation is carried out with these more accurate boundary values 


and so on. This iterative process converges rapidly. 

5. Compressible flow 
In this section incompressible values will be distinguished by a suffix ° 
We have the following equations 


(pur), per) _ 4 a 


Continuity: 
Cx cor 
. 7 ‘ : dp AY) 
Bernoulli: q dq a*dp/p. ie. dL = —_. (2! 
M?*p 


Vorticity: 
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Writing w qcos@ and v qgsin @, this equation becomes 
ou Cc Cc ° oL . 
r|—sin@ cos 0 cos 6 sin 6 —]-+-siné 
C2 Ca cr cr 
oL _ oie ; 
rM2( cos 6- 1 sin@ } rX’. say. 22) 
Cx cr 
\lso (21) ean be written 
{ o6 é : of oL rc ' 
r {cos 6 sin 6 sin @ cos 6 — >. (23) 
Ca Ca cr cr , q 
. {99 »)2)\ ; 
rom (22) and (23) we find 
( oL. 0 ; a. 
r| cos 6sin@ rid’ cos @4 sin 8), 24) 
ca co? qd 
. coo a — s C ~ 
ana | sin-é r{A’ sin @ “cos 4). (25) 
Ca or q 
We now transform to the (¢,%) plane using 
4 € c 
u . and v; 1 1g mn, 
? Ch Ou cr Ch Cus 
which follow from (1), then after some calculation (24) and (25) yield 
( L coo « ° ’ 
; : A, COS x sin x, (26) 
C@O Ww qq; 
ale) cL : C ie 
7 A, sina COS x. (27) 
Ch C2 qq, 
: oL : cL sind 
where A, M*\ cos ' SIN «a (28) 
\ CO ous rq; 
Cross differentiation of (26) and (27) yields 
¢ OL LoL (1  — | 
1 . | . [A, cos y+ — SIN \} 
Cw Cw @\? Ch | } 149; | 
C ° 
- I), sin .——*-cosal, (29) 
Cus | 99; 
o {1 o6 of\ { | ; 4 
—\- q (7 a A, SIN a — COS ) 
eb\ r ad ab i 99; | 
C : 
/A, cos x4 sin ey (30) 
cas \ 19; | 














364 L. C. WOODS 


which correspond to equations (24) and (25) of reference 1 for two-dimen. 





sional motion. then 
If we ignore «, the error introduced will be small for subsonic flow, ang I The | 
further, if the vorticity is zero, (29) becomes but i 
0 [( eL 0 /leL 0 (M2 oL é [sin @; oy, [tweel 

“al =i) : Zale 28) ‘al r =) taal =} a, 


If (8) and (10) have been solved, r, q;, and sin @; will be known as fune. _E 
tions of ¢ and y%, and hence (31) can be solved by relaxation. 


ae ‘ : ; soe te Th 
A finite difference equation representing (31) is easily found. Referring 
; ‘ genel 
to Fig. 1, we find that The | 
e 
] ] 
r,(L,—L,)—1,( Ly>— L4) + — (L,— Ly) ——(Lp— Ls) press 
Ya "h = 
avoic 
M? M? — (sin 8 sin @ a 
(=) ( L, Ly)— : ( Ly Ls) oe > ie al ’ a 
r Ja Y Jo Gl sa yee the a 
and thus Th 
1— UM? 1— M? theor 
=) beni bak. 2 
r d - rh dictes 
1M) (1—M? _(siné\ (sind) —_,, pont 
latte Tt sl i Ly ce eae —s]> \4I Band t 
| r d r b arya  \UT Jo 
; ss ; . wall | 
is a suitable difference equation representing (31). that 
The special case r = 0 requires some attention. On r = 0, re 
ulvel 
é L -0. 8 0 r C (7, @) op “ 0. will n 
Cus Cus (dh. bs) Cx rq; tis 
ye : ‘inite 
and oe ~—— a “i 0. * t 
Ch O(d.us) Cx qi In 
_— ’ there 
Substituting these special values in the expanded form of (31) we find mn 
aL c M2 . L (33) irst « 
cd? od ch n ac 
— the p 
6. The boundary conditions ‘hk. 
nuit 


If R is the radius of curvature of the profile in the (a, 7)-plane. and ns 
the direction normal to the profile boundary, then 


oi q : 
on R’ 

Hence, ob a , ! é 
on R 


and since 
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eL | 
cp Rq,r 


The boundary conditions are, of course, not specified for the (¢,%)-plane, 


but if the incompressible flow has been determined the relationship be- 


tween the (x.7) and (¢,y) planes will be known, and hence the boundary 


nditions can be immediately deduced for the (¢,) plane. 


7, Example: the design of a contraction 

The problem of this section arose in connexion with the design of a 
seneral-purpose wind tunnel which may shortly be erected in New Zealand. 
[he problem was the design of a contraction ratio of 6-83 with a linear 
pressure gradient along the wall. Adverse pressure gradients were to be 
voided if possible. The ratio of contraction length to the radius at the 
rgest end was to be about 3-6. Air velocities were small enough to permit 
the assumption that flow was incompressible. 

The solution is shown in Fig. 2. The curve A BC is that deduced by the 
theory of sections 3 and 4. It deviates considerably from the curve pre- 


icted by the linear theory. At first the contraction was started at the 


ipoint B, where ¢d = 0, but to give a more gentle lead into the contraction 


nd to avoid too serious an adverse pressure gradient on the straight tunnel 


wall preceding the point B, the wall was ‘relaxed’ back to the point A so 


hat constant pressure was maintained in the section AB. A very slight 
dverse pressure gradient still survives in the section AD, and in fact this 
vill not disappear until the contraction is of infinite length. In other words, 
t is not possible to avoid an adverse pressure gradient in a contraction of 
nite length. This becomes evident when equation (14) is applied to a 
int such as F in Fig. 2. It is clear that AG > EF, i.e. from (14) q4 < qp, 
herefore the pressure at A is greater than the pressure at F’, and so there 
san adverse pressure gradient along AD. Mr. W. Rainbird of Cranfield 
rst drew the author’s attention to this point. There must similarly be 


n adverse pressure gradient in the section CF. It will be noticed that 


he pressure curve has been rounded off at B and C to avoid a discon- 


nuity in the gradient. Table 1 gives the coordinates of the contraction. 


TABLE |] 


Coordinate S of the Contraction 


‘96 | III5 | O15 g'17 8°56 8-14 
‘ 


406 | 44°6 
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The calculations were made as follows 


1. An approximate solution was found from the linear theory. 





If the velocity at the inlet to the contraction is unity then, to give a 
2 | linear pressure gradient, the velocity at a point x, in the contraction of 
enoth XY. must be given by 
| r 


- J q° | + —(k*—1), 
X 
vhere k is the contraction ratio, and the origin is taken at the entrance to 
cal the contraction. Now 
| xX X : 
F i ’ 1 2 & x 5 
=F d vy da if — (e* 1)\2 dx 2 1+ ~ (k?—1) Pn ie ° 
Sad ; } X 3 k2@—1] X 
0 0 (35) 
When x X ri) 3X (k?— 1). 
no | 


Suppose that at 2 = 0,r = %, then X = 3-6y (specified ratio X /7~ = 3-6), i.e. 


( 


d = 23-6 x (6-832—1)4, = 16-27. 


= Now % on the boundary is given by 
4 
® To 
Wy, rar = Lre. 
a é . 
| sa Hence, ds, / %/32°4; e.g. if % 16. us,/h 4. 
| ha ial rr ; 
| =] Suppose that there are 10 divisions of ¢ in the contraction, then there must 
oO = 
—~ be 5 divisions of ys. These divisions are shown in Fig. 2, and the above argu- 
|, | ment shows why 16 was selected as the value of 1. 
——* ‘ : : , ' 
‘egardless of the nature of the pressure gradient, the linear theory yields 
3) 
—— 2 2 
| ee {o"0> 
| . . 1 92 
lv Le. r, 1/935 (36) 
4 ind it was from this equation that the dotted line in Fig. 2 was calculated. 
2. Using (35) we have | 27 
ra) d(x), (37 


——) and hence from the given pressure distribution. p p(x), we can deduce 
a relationship q? q°(h) to use in equation (16). Further, from (36) and 
— or 


) 


) we find approximate boundary values of f 
f f(d) f{{dh(a)! y® == sila. 
Equation (13) shows that we merely need to fill in the values of f along each 


os _| equipotential line so that the gradient is linear, with f = 0 on 0, and 


f(¢) on & = y,, to obtain a fair approximation to the f field. 


50 


3. The residuals are now calculated and the relaxation carried out. A 
pomt to notice here is that if the equipotentials are correctly spaced for a 


reasonably small mesh size (to avoid large errors in the neglected terms of 
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the exact difference equation) near the small end of the contraction, they 
will be far too largely spaced at the large end (see Fig. 2). This can be 
easily remedied by inserting extra equipotentials at the large end so that 


The mesh actually used is shown in Fig. 2. 
. Next equation (18) was used to integrate along the boundary to yield 


x(¢). 


5. Steps 2 to 4 were repeated twice, after which no significant change 


urred in the boundary conditions, and the f field was thus ‘settled’, 


6. Difference equations representing (9) and (10) were then used on the 
f field to complete the solution. 
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Design for minimum length of contraction ratio 
It is of some interest to be able to design a contraction so that, given the 


contraction ratio and the maximum permissible adverse pressure gradient 
(c), the contraction is of minimum length 1. Referring to Fig. 3 we have to 
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determine the pressure curve AB so that, given e,/ isa minimum. This is 
» difficult problem, but it could be solved by relaxation. Consider the linear 
sradient family shown in Fig. 3. The experience gained with the problem of 
section 7 leads the author to believe that of this family the one labelled 
{CDB (CD extends over one mesh interval only) would result in the 
shortest value of /. Theoretically the gradient of CD could be made infinite 
mit, as in the corresponding two-dimensional problem, this would result in 
in infinite spiral at C'D (@ becomes infinite at a discontinuity of log(1/q)) 
which would be impracticable. If it is true that ACDB (with CD vertical) 
s the optimum curve of the linear gradient family, then it appears that 
the significant point for minimum length is that the constant gradient 
portions, AC’ and DB, should be as long as possible. It also appears that 
ACDB is the optimum curve for all possible families. A series of calculations 
by the method of this paper could prove or disprove this point. If ACDB 
proved to be the optimum curve, then it would be necessary to carry out 
. further series of calculations to establish the optimum value of x/l at which 
CD is to be placed. At values of x// other than this optimum it is likely 
that the adverse pressure gradient at one end of the contraction would not 


reach the permissible maximum e. 


APPENDIX 


Since this paper was written a simplification of the relaxation treatment of equa- 
8) has occurred to the author. This equation can be written 
*log( o? 
J) FF 6. (38) 
od” Cus? 
with the same order of error term as for the difference equation (11) we can 
38) (see Fig. 1) in the forn 
| } log( J,) log if 2 log( fo)} i Jo iP 2h} Xo 
Mt _ . va , ‘ 
\—) log( fi fs/f2)+fet+hi—2fe = Xo- (39) 
n 
ithe n coefficients | 
\ oX m)\* oX { m)\? /. ) 
1, = =f has = 2\1 ae (40) 
ahaa (nd [iw = PUG) | 
residuals are now calculated from (39) and relaxed using equations (40). Equa- 
s (39) and (40) have obvious advantages over (11) and (12). They are simpler 
avoid the use of the intermediate points b and d. 
With this m dification, (16) and (18) become 
m\* (4m2 (m)\? (f,—fe)?\4 
X (—) log( fi fs/f3) +2 fi—fo) +2 | (—j <2 Is (41) 
} q n to 
} 1/m\? , oo 
t | ly Ja) -=\ ) log(hi fa/fo) |» (42) 
} o vt 


spectively 








A NEW RELAXATION TREATMENT OF FLOW 





REFERENCES 
1. L.C. Woops, The Design of an Aerofoil Profile to give a Specified Velocity Distriby- 
tion at a Given Stream Mach Number, A.R.C. F.M. 1412 (Nov. 1949), 
2; A New Relaxation Treatment of the Compressible Two-dimensional Flow about 
an Aerofoil with Circulation, A.R.C. 13,034—F.M. 1429 (March 1950). 
. R. V. SoutHWELL, Relaxation Methods in Theoretical Physics (Oxford, 1946), 
4. A. Tuom, Arithmetical Solution of Problems in Steady Viscous Flow. R. & M. 1475 


(May 1932). 
G. K. BATCHELOR and F. 8S. SHaw, A Consideration of the Design of Wind Tunnel 


Contractions. Report ACA. 4 (March 1944). 


5. 














\PPLICATION OF RELAXATION TO THE ROTATIONAL 
= FIELD OF FLOW BEHIND A BOW SHOCK WAVE 


By A. R. MITCHELL (Vathematics Department, United College, 


a). The University, St. Andrews) 
MI. 147 
Received 6 February 1951] 
1 
SUMMARY 
\ development of the relaxation technique due to R. V. Southwell for evaluating 
ibsonic isentropic regions of compressible flow in two dimensions is used to evaluate 
nixed subs c-supersonic regions where the flow is rotational. 
The new technique is employed to determine the complete field downstream of the 
bow shock wave formed when a parallel supersonic flow of Mach number 1-8 impinges 
na blunt-nosed two-dimensional obstacle. The vorticity is caleulated everywhere 
n the field and it is found that the stream-line starting at the point of maximum 
rticity on the bow shock wave is a locus of high vorticity in the rotational field 
lownstream of the shock wave. 


| 1, Introduction 

Two attempts have been made so far to obtain a numerical solution of the 

field of flow behind a shock wave. Emmons (1) developed a method of 

relaxation, applicable to subsonic regions only, which he used to evaluate 

the rotational subsonic field of flow behind an arbitrarily placed shock in 
» hyperbolic channel. Maccoll and Codd (2) investigated the important 





practical problem of a supersonic flow of Mach number 1-5 streaming past 
a flat-nosed obstacle, the position of the bow shock wave having been 
obtained from photographs. They examined the subsonic region behind 
the shock wave by a relaxation method, the solution being based on the 
assumptions that the flow was irrotational and that the sonic line started 
at the obstacle corner and intersected the stream-lines at right angles. 

In the present paper Southwell’s relaxation method (3) for evaluating 
subsonic irrotational regions in two dimensions, as extended to cover 
supersonic regions by Rutherford and the present author (4), is further 
extended to evaluate mixed subsonic-supersonic regions where the flow is 
rotational. The new technique is then used to determine the entire field 
downstream of the bow shock wave formed when a uniform parallel flow 
of Mach number 1-8 impinges on a square-nosed two-dimensional obstacle. 

Several attempts have been made to obtain an analytical solution in the 
hodograph plane of the two-dimensional field of flow behind a bow shock 
wave. These investigations, reviewed by Guderley (5) and Busemann (6) 
all deal with free stream supersonic Mach numbers close to M 1. This 
restriction on Mach number permits the assumption of approximately 


[Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 3 (1951)] 
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isentropic conditions behind the shock, and so Chaplygin’s form of the 
compressible flow differential equation in the hodograph plane is used, 
In a very recent paper Goldstein and Lighthill (7) gave some indication 
of the difficulties likely to be encountered in using the hodograph plane to 
obtain analytical solutions of flows with vorticity present. In the simple 
example of the flow of an incompressible stream with uniform shear past 
a circular cylinder, the hodograph plane was shown to be a Riemann 


surface of six sheets. 


2. The fundamental equations 
For the two-dimensional steady motion of a rotational non-viscous 


compressible fluid with uniform conditions at infinity, the stream-function 


us must satisfy (8) 
O [1 eb O {1 ed » oS 
C oy i Cus , } ; _— (1) 
Cx\p ox} ' cy\p cy)  B eds 


In this formula p, p, and S denote the density, pressure, and entropy 
respectively, while # is the gas constant. After some reduction this 
equation may be written in the form 
V2(xp)—pV2x 4+ 7 at 2) 
, Ry by 
where x = p 
In the flow behind a shock, the conditions, although adiabatic, are not 
necessarily isentropic. However, the entropy maintains a constant value 
along each stream-line, and so Bernoulli's equation for the rth stream-line 


downstream of the shock may be written 


(x). .. l, P (ps)r\ 
cl yl op (Dodd 
where wu is the speed, y denotes the ratio of the specific heats, and c, is the 
speed of sound at a point of stagnation; also (p,), and (p,), are the values 
of the stagnation density and pressure respectively on the rth stream-line 
downstream of the shock. The pressure may be eliminated by using the 
adiabatic gas law 
Sette & a ON ‘ 
pY — pk (p.)} 
where Sp, p,, and p, are the entropy, stagnation pressure, and density on 
all stream-lines upstream of the shock, and S, and k, are the entropy and 
adiabatic gas constant on the rth stream-line downstream of the shock. 
Bernoulli’s equation thus becomes 


(“) = | ‘ AS] } 
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which may be written in the form 


Z iF | x ¥ ? | Le - (~\) (4) 
1 | Xs) | ce \\ ex a) " 
(2) and (4) are the fundamental equations for the flow behind a shock wave. 


It will be advantageous to write these equations in the non-dimensional 


forms p l oS 


V>( Ks) pV7x 0 (5) 


Bx yve Cus 


2 oa f(b? (eb?) 
‘ d f] 2(1—} , 2 al(: | ( (6 
wn 7 x VX \\ a5 ay | >) 


To obtain these equations, some significant linear dimension / pertaining 
to the problem under consideration is chosen, and the flow parameter on 


the rth stream-line behind the shock is expressed as 
G 


(Ps)p Cs h’ 


in which G is the mass flow per second under free stream conditions. With 
. slight change of notation, x, ws, p, x, y now stand for the non-dimensional 
quantities y/(x,),.u/G@, p/(p,),,«/h, y/h. Eliminating the pressure by means 


ot (3). (5) he comes 


Ta ; l Los 
V2(xu) pVex T > (ny l , 0. ( 
ye x27) B ous 


~I 
— 


Thus the fundamental equations (2) and (4) take the non-dimensional 
forms (6) and (7). 

Before the last two equations can be used to evaluate the flow, it is 
necessary to obtain the values of v, and (é@S/as)/A on every stream-line 
downstream of the shock. As a consequence of the fact that G/h and c¢, 
have the same values on both sides of the shock, the flow parameter ratio 


for the rth stream-line is given by 


ve os 

v= (pe 

where v is the flow parameter on all stream-lines upstream of the shock. 
Now using the gas laws 


2 — YPs_ Y\Psh Pe_ pp _ Pr 

; p (Ps), px pi 
and \Ps)r k. = Po 
(P.)p Ps 


where the suffixes 1 and 2 indicate quantities at the shock on the upstream 
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: : at 
and downstream sides of the rth stream-line respectively, the flow para- 
a | 
meter ratio becomes 
erl 


Ve (22) ; (?)” = | 
ve Np Ps 
or in terms of the ratio p,/p, only, 
Ve (Pa)? (y—V)(Po/ p+ (y+) Pre 
(7: bs + 1)(p2/Py)+(y 4 


It will be shown in section 4 how (8) together with the standard shock- 


, (8) 
vy” 


wave data can be used to evaluate v. and (éS/éd%s)/A on stream-lines 


’ 


downstream of the shock in the problem of this paper. 


3. The relaxation method 
In order to use the method of relaxation, the field to be examined is 


covered with a regular network. In the problem considered in this paper 
an 
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Fic. | 
a square network is used. The relaxation method and the use of patterns wi 
to eliminate residuals are discussed fully in (3) and need only be sum- 
marized here. In all calculations the value 1-4 has been taken for y. TI 
The finite difference approximations to (6) and (7) for a square net of on 
non-dimensional mesh size a are, in the manner of (3), an 
R —4  .—aay) = (Mra | ye 1 uf b.— 2b)! (9 In 
0 Xo Xo ‘ (y ) at | "OT & (ys ;,— 2 sg) | 7 ne 
i=1 
Tes 
2 " 
_ , as 1 {eS . 
and KF, > x(b;—) ; 0, (10 of 
a? 2)  ~@yt1) g 
— V(r )0 Xo” B\ cf} 
where the suffixes 0, 1, 2, 3, 4 are as shown (Fig. 1). Thus for a given of 
W-distribution, using (9), Ry and hence y, can be calculated at every node ol - 


the network. Equation (10) will then enable the residual F, to be evaluated 
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at every node of the field. Modification of the original £-distribution with 
a view to eliminating these residuals is carried out either by trial and 


error or by a relaxation pattern. 








In the present case the changes in bF; (j 0, 1, 2, 3, 4) consequent on 
A 
4 Y, 
rA 
D Yr / 
B O : D vy; Wn 








Fic. 2 


an assumed modification 6d, of d& 


" 4, are given by the relaxation pattern 


formulae 


oF. SY 0 ‘ a f’(R)i (bo—,)?4 
font za 
(2y+1)(y—1) 1 (eS) 5g, . ; 
y x4 } 5 ) c fy (R)}o( 4yho - > ¥;) dbo 
XT ang BV Ab), io 
(11) 
und 
; (y—1)(v?), , . 
SF | Yo Al f’( R) (bo b;)( dbo - > w;)— 
wa" i=1 
(2y L)(y 1) 1 foes ie R) 1 »— 1,2,3,4 
2271 Blea) ( )} {(Yo—;) Of (J * Enact 
‘ , (12) 
where, as in (4), f'(k x 


A{(yt+1)x7%1-1/2—1} 

These formulae are also obtained in the manner of (3). They are, however, 
only approximately true, as it was assumed in their derivation that (€S/0xs), 
and (v2), did not vary with the increment 8, of %. The consequent 
inaccuracy in the relaxation pattern formulae is unimportant except in the 
neighbourhood of the speed of sound, where a trial-and-error method of 
residual elimination must be used on another account, detailed in section 4 
of the present paper. 

Because of the curved nature of the bow shock wave, which is a boundary 
of the region examined in the problem of this paper, the field contains 
several irregular stars. The evaluation of such nodes is discussed in (3). 
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The method found to be most suitable in the present problem is illustrated, 
O is an irregular star (Fig. 2(a)). The four associated points are A, B, (. 
and D, where OB < OA < OC OD =a, the non-dimensional mesh 
size. A’, C’, and D’ are such that OA’ OC’ OD’ OB = b, say, 
where } is non-dimensional. The stream-function &% is known at A, B. (. 
D, and O, and so by interpolation the value of % is found at A’, C’, and D’. 
Thus the irregular star centred on O has been replaced by a regular node 
of a square net of non-dimensional size } (Fig. 2 (b)). 


4. The problem to be solved 

The particular problem described in this paper is that of a uniform 
parallel stream of Mach number M@ 1-8 (v 0-40) flowing past a square- 
nosed two-dimensional obstacle. The shape of the bow shock wave was 
obtained from a photograph made in the N.P.L. (9), and the extent of the 
field considered is from the axis of symmetry to the first stream-line which 
is not appreciably deflected by the shock. At such a stream-line the shock 
angle should be approximately the Mach angle sin-!(1/1-8). It is required 
to find the complete pattern downstream of this shock. 

With a view to applying the methods of the preceding section, the 
significant linear dimension h, which was introduced in order to obtain 
non-dimensional quantities, was chosen to be five semi-widths of the 
square-nosed obstacle. The first stream-line not appreciably deflected by 
the shock is at a distance 10f from the axis of symmetry. The shock 
cuts this stream-line at A and the axis of symmetry at C. The sonic 
velocity point on the downstream side of the shock is B, the stagnation 
point is D, the obstacle corner is LH, and the perpendicular from A meets 
the obstacle side at F. The methods of the previous section are now used 
to evaluate the field inside ABCDEF, where the rotational flow is mixed 
subsonic and supersonic. The problem is illustrated (Fig. 3). 

The boundary conditions for this region consist of a knowledge of the 
stream-function y%, the Mach number , and the slope of the downstream 
stream-lines along AC, and also of % along CDEF. It is also necessary at 
this stage to obtain v? and (éS/és)/A at the shock on the downstream 
side. At every point of the shock wave the tangent makes the shock angle 
with the direction of the axis of symmetry. It is well known that for a 
given Mach number and deflexion angle there are two possible shock-wave 
angles. The larger angle corresponds to the strong shock and the smaller 
angle to the weak shock. For example, at points C and A on the bow 
shock wave where the angle of deflexion is zero, the shock conditions are 
strong and weak respectively. If G@ is the point on the bow shock wave 
where the stream-line experiences its maximum deflexion, then all points 
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ated from C' to G correspond to strong shock conditions and all points from @ 

B, ( to A to weak shock conditions. There is, of course, only one shock-wave angle 
mes] possible for the maximum angle of deflexion, since the weak and strong 
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smaller shocks coincide at G. From the oblique shock chart of (10), the pressure- 
1e bow ratio p, p, can be read off for an initial Mach number of 1-8. Substitution 
yns art of this pressure ratio in (8) will then give the flow parameter ratio at every 
k wave point of the shock wave. The entropy has a constant value S, everywhere 


points in the field upstream of the shock, but takes a different value S,. on every 
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stream-line downstream. The entropy increase (S,—S,) is obtained at all 


points of the shock from (11). Since the stream-function % is also known 


at every point of the shock, we may plot (S,—S,)/7 against ys and read off 
the required quantity (@S/es)/7 at all points on the shock. v?/y? and 


r 


(CS/0s)/A are each plotted against ib (Fig. 4). 























The curves v2/v? = constant and (éS/és)/ 2 = constant are the stream- 
lines, and so it follows that Fig. 4 may be used not only for points on the 
shock wave, but also for all points downstream of the shock in the field 
to be examined. 

As already mentioned, the field inside ABCDEF was covered with a 
square net of non-dimensional mesh size a = j}. Now it was shown in (4) 
that the only difficulties encountered in employing the relaxation method 
to mixed subsonic-supersonic problems occur in the neighbourhood of the 
sonic line. The method outlined there for overcoming such difficulties was 


employed in the present problem in the following way. 
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with shape of the sonic line is known, BCDE will constitute a subsonic region 
. in (4 and ABE F a supersonic region. Using (9), R, can be calculated from the 
catia v-distribution at every node of the field. From Fig. 1 of (4) it is seen that 
of the ior each value of R, there are two possible values of x, a value less than 
an ine |:25 if the node is in the subsonic region, and one greater than 1-25 if the 


node is in the supersonic region. Although the sonic line is not known 
initially, the appropriate ¥» can always be selected at each node. Thus F, 
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is calculated at every node using (10). For the nodes near CD with com- 


paratively low subsonic velocity and for the nodes near A F with compara- 
tively high supersonic velocity, the modification in % is made solely with 
a view to the annihilation of the residuals F,. At such nodes the relaxation 
pattern formulae (11) and (12) may be used in the elimination of the 
residuals. However, on approaching the sonic line from either side it is 
found that a modification of y% often results in a value of R which is greater 
than the maximum possible one of 0-068. In such a case there is no 
possible value of y. Consequently in modifying the value of % at nodes in 
the vicinity of the line BE with a view to eliminating the residuals, the 
value of R must be kept below 0-068. At such nodes the pattern formulae 
do not yield a definite result, and a method of trial and error is adopted. 
The sonic line is eventually obtained in the manner of (4). The complete 
field of flow showing the lines of constant Mach number is illustrated 
(Figs. 5 and 6). 

In the flow of a subsonic compressible stream past a wedge (4) it was 
found that the sonic line met the wedge just in front of the apex. Conse- 
quently the field close to the obstacle corner in the present problem was 
examined in order to determine as accurately as possible where the sonic 
line met the obstacle. The net required to produce such information had 
a non-dimensional mesh size, a = ;},. It was found that the sonic line 
again met the obstacle just in front of the corner. As a detailed examina- 
tion of the corner is not the primary object of this paper, the mesh size 
was not reduced further. 


5. The rotational field 
At any point in the rotational field of flow behind the bow shock wave, 
the vorticity is given by (8), 


p os 
W . 
BR ols 
This equation can be written in the non-dimensional form 
P l os 9 
@= — (13) 
yv, B eds 


where with a slight change of notation, w, p, and % now stand for the 
non-dimensional quantities (h/c,)w, p/(p,),, and %/G respectively. 


s 


Now, as mentioned earlier, v?/v? and (@S/és)/#A depend only on y and 
so the quantity (@S/és)/yv,Z has a constant value along every stream- 
line behind the shock. The non-dimensional quantity p, like the Mach 
number WM, is a function of y at the point in question, and so (13) enables 
the non-dimensional vorticity to be calculated everywhere in the field 
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downstream of the shock. Lines of constant vorticity are shown in Figs. 5 
and 6, along with the lines of constant Mach number. 

On the shock wave itself, the value of the vorticity does not depend on 
the field downstream of the shock. Fig. 4, constructed from oblique shock- 
wave data, enables w to be calculated from (13) without reference to the 
field downstream of the shock. Because p/v, decreases gradually with 
increasing ys along the shock, the vorticity is found to be a maximum just 
before the point of inflexion on the graph of (S.—S,) against %, which was 
obtained from (11). This point on the shock wave where the vorticity is 
a maximum is of some importance, since the stream-line starting there will 
always pass through regions of comparatively high vorticity. 

A formula will now be derived for the angle between the stream-line and 
the constant Mach number line through any point in a rotational compres- 
sible field. Consider the field covered by two orthogonal families of curves 
|= constant and n = constant, the latter family being the stream-lines. 
The velocity q at a point is given by g = q(/,n). Thus along a line of 


constant velocity. ] 
(‘| (4) ( n 0 
ol , Cn}; dl q ; 


If d is the angle between the constant velocity line and the stream-line 
through a point in the field, it follows that 


tan [> dn - h, (cq él), 
=m h, il), h, (éq/én),’ 


where h,d/ and h,, dn are the elementary arc lengths along the stream-line 
and its orthogonal trajectory at the point. Now the vorticity at a point 


is given by l (eq q 
7 rat i) L’ 


where L is the radius of curvature of the stream-line at the point. Elimina- 
tion of (ég/én),/h,, leads to the result 
L l (4) 
q cs n 
Lw ; 
q 


ls is the elementary are length taken along a stream-line. 
Since ¢ 


tand 


where ¢ 
has the same value on every stream-line behind a bow shock, 
it can be deduced from Bernoulli’s equation that a line of constant velocity 


is also a line of constant Mach number. In addition, it is easily shown that 


l (ec | /oM 
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The angle between the line of constant Mach number and the stream-line 


‘oM 
| cs n 


z (1+ ar) ; o 147 ar} 
L 2 2 


where all the quantities in (14) are non-dimensional. On the obstacle sides 


through a point is therefore given by 





tand 


in the present problem, the radius of curvature is infinite and the vorticity 
zero. Thus the constant Mach number lines meet the obstacle everywhere 
at right angles. 

[t is worth while pointing out that in the rotational field of flow behind 
a bow shock, since c, is constant everywhere, a line of constant dimensional 
velocity is a line of constant Mach number, constant non-dimensional 
pressure, and constant non-dimensional density. However, since (p,), and 
(p,), Vary with r, lines of constant dimensional velocity are not lines of 


constant dimensional pressure or constant dimensional density. 


6. Conclusion 

The most significant feature of this paper is the adaptation of the relaxa 
tion technique to the evaluation of mixed subsonic-supersonic fields 
changes in entropy being taken into account. The only difficulties ocew 
in the neighbourhood of the sonic line, and these are overcome by the 
method outlined in (4) and in section 4 of the present paper. 

The importance of taking into account the rotational nature of the flow 
behind the shock wave in the problem of the present paper is illustrated 
in Fig. 5, where it is seen that far downstream of the shock, where the 
flow is approximately parallel, the Mach number varies from 1-42 on th 
obstacle side to 1-80 on the first stream-line which passes through the shock 
without appreciable deflexion. If the flow behind the shock had beet 
assumed irrotational, the Mach number of the flow far downstream would 
have been approximately 1-80 on every stream-line. Thus in the present 
problem the assumption of isentropic conditions behind the shock would 
have given a very inadequate approximation to the true field of flow. Evel 
in the problem discussed by Maccoll and Codd (2), where the free streat 
Mach number is 1-50, the flow cannot be considered even approximately 
irrotational behind the shock, as the gain in entropy at the shock i 
sufficient to produce a parallel flow of Mach number 1-36 on the obstacl 
side far downstream of the shock. 

The point on the shock wave at which the vorticity is a maximul 


is of major importance in relation to the rotational field of flow down: 
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stream of the shock, since the stream-line starting there is a locus of high 
vorticity. 

Finally, the formula derived in this paper shows that the angle between 
the stream-line and the constant Mach number line through a point in a 
field of rotational flow is a right angle only if the stream-line has no 
curvature and the vorticity is zero at the point. This is the case on the 


obstacle sides in the problem of this paper. 
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\ NOTE ON KLITCHIEFF’'S PAPER ON “BUCKLING 
OF CONTINUOUS BEAMS ON ELASTIC SUPPORTS’ 


By A. CHARNES (Carnegie Institute of Technology, 
Pittsburgh 13, Penn., U.S.A.) 


[Received 13 November 1950] 


IN a recent paper (1) J. M. Klitchieff arrives at an equation (5’) for a 
function Q which may be written 
mQ > {[(2k+1)m—rP—m?}-{(2k+ 1)m—r}-* 
k= 

He is interested in the function ® = m'z-*Q-! for r = 1 and m 
and by approximating the sum of the series obtains the formula 

@ = 27-2m4(5m?2+ 1)(m2—1)-2(4m?2—1)-1. 
Professor E. Saibel has drawn my attention to the fact that he obtained 
by a different method the following expression for ®, viz. 

® = 21+ cos(z/m)]-?. (2’) 

The two results contained in (1) and (2’) can be shown to agree. For the 


series (1) can be summed by residues (2) to give 


® = 4-1 cosec2a( | 4 | cot a a. nid —cot = = » (3) 
2m 27 2m 


which for r 1 becomes 


= 4 cose? 7% "), (4) 


2m 


and it is easily shown that (2’) and (4) agree. 
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